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A  BSTRA  CT 


This  paper  deala  with  the  formation  of  eddiea  in  a  straight 
parallel  or  zonal  flow  and  with  the  .subsequent  modification  of  the 
flow  profile.  The  fluid  is  taken  to  be  homogeneous  and  inviscid. 
Numerical  analogues  for  the  physical  equations  are  developed  in 
detail  and  are  analyzed. 

The  work  begins  with  the  linear  theory  of  dynamic  stability, 
Ntanerical  analogues  are  developed  to  determine  the  evolutioxi  of 
perturbations,  sinusoidal  along  the  flow,  which  are  initially  pre¬ 
scribed  with  arbitrary  wave  number,  amplitude,  and  tilt  varia¬ 
tions,  and  which  are  superimposed  on  arbitrary  flo>*'3.  These 
flows  are  straight-parallel  and  are  unbounded,  or  are  lialf-bounded 
or  bounded  by  plane  surfaces.  Integrations  are  carried  out  for  an 
unbounded  flow  profile  with  an  inflection  point.  Unstable  perturba¬ 
tions  are  isolated  and  the  unstable  spectrum  is  determined. 

A  numerical  analogue  for  the  finite-amplitude  problem,  by 
which  one  can  study  the  transfer  of  energy  from  the  mean  flow  to 
the  eddy  then  developed.  The  most  unstable  perturbation,  linearly 
determined,  is  taken  as  a  small  but  finite  disturbance.  The  integra¬ 
tion  is  carried  out  and  reveals  the  continued  growth  of  the  eddy  and 
the  modification  of  the  mean  flow. 

This  method  of  investigation  with  added  lapse  rate  and  com¬ 
pressibility  is  discussed  as  an  approach  to  turbulence,  and  to  the 
modification  of  wind  shear  and  lapse  rate  by  the  developed  eddies. 
The  general  problem  of  numerical  analogues  fcr  integrations  re¬ 
quiring  finite  time-steps  is  also  briefly  discussed. 
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A  NUMERICAL  INVESTIGATION  OF  THE  BAROTROPIC 


DEVELOPMENT  OF  EDDIES 


1.  INTRODUCTION 

The  study  of  the  beha\'ior  of  sniall  disturbances  superimposed  on 
straight  parallel  flow  originated  in  classical  hydrodynamics  in  connec¬ 
tion  with  the  stability  of  laminar  flow*,  {For  a  review  of  the  literature 
and  an  excellent  bibliography  the  reader  is  referred  to  Lin.  )  In  the 
present  investigation,  we  treat  two-dimensional  motions  of  an  inviscid 
homogeneous  Uuid. 

A  straight  parallel  flow  is  said  to  be  dynamically  stable  if  all  super¬ 
imposed  infinitesimal  perturbations  aie  out  or  remain  bounded  at  all 

tirnea.  Otherwise  the  flow  is  unstable  and  the  profile  must  be  considered 

8 

a  transitory  state.  Rayleigh  showed  that  an  inflection  point  in  the  flow 
profile  is  a  necessary  condition  for  instability.  Sufficient  conditions  have 
been  established  for  flow  profiles  of  certain  general  types;  however,  no 
complete  theory  has  yet  been  developed. 

9 

Considering  flow  profiles  with  a  point  of  inflection,  Tollmten  showed 
that  fc  ‘  symmetrical  profiles  the  existence  of  a  neutral  oscillation  im¬ 
plies  a  transition  from  stable  to  unstable  solutions.  In  an  investigation 
of  unbounded  broken  profiles,  Holmboe^  iound  that  those  having  maximum 
vorticlly  at  the  inflection  point  are  unstable,  whereas  those  having  rnini- 
mum  vorticity  at  the  inflection  point  are  stable.  This  condition  is  appar¬ 
ently  not  sufficient  for  flow  in  a  channel.  An  example  which  does  not  con¬ 
form  is  a  sinusoidal  profile  in  a  channel  whose  width  is  less  than  half  a 
wave  length  of  the  profile.  In  this  system  the  boundaries  inhibit  the  de¬ 
velopment  of  unstable  disturbances.  However,  if  the  channel  is  wider 
than  half  a  wave  length  the  current  is  unstable. 

(Author's  manuscript  approved  ?.4  October  1958) 
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In  the  present  work  a  method  of  testing  the  dynamic  stability  of  arbi¬ 
trary  profiles  is  developed.  Initial  value  problems  are  solved  by  direct 
numerical  integration. 

2.  THE  SYSTEM 

The  fluid  is  homogeneous,  incompressible,  and  nonviscous.  The 
model  is  two-dimensional  in  the  sense  that  there  is  no  motion  -and  no 
variation  in  any  of  the  fields  or  boundaries — along  the  third  dimension. 
The  effect  of  gravity  in  the  model  is  trivial  since  we  sliall  not  be  dealing 
with  any  free  surfaces,  and  the  gravity  field  is  herewith  dismissed. 
However,  for  orientation  we  may  refer  to  the  straight  parallel  flow  as 
being  horizontal,  parallel  to  the  x-axis,  with  the  speed  varying  in  the 
vertical  along  the  z-axis,  although  some  applications  of  the  models  to 
meteorological  problems  may  be  otherwise  orientated. 

The  total  velocity  field,  ,  is  considered  to  be  the  resultant  of 

two  component  fields.  One  of  these  is  the  straight  parallel  flow,  t  Jq.  , 
which  would  be  a  steady  state  if  it  existed  alone.  The  other  component, 
V”  ,  shall  be  called  the  disturbance  or  perturbation  velocity.  Only  in 
special  cases  is  the  resuHant  field  a  steady  state.  In  general  the  flow 
evolves.  Its  evolution  shall  be  absorbed  by  "V,  the  straight  parallel 
flow  being  maintained  constant  by  choice.  Thus 


'V  =  U(2):  -H  V  CX.H.t) 


(1) 


Because  V  is  the  difference  between  two  nondivergent  fields,  we 
may  represent  it  by  a  stream  function. 


v=  jxVV-ua-h  urik. 


(2) 
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The  total  vorticity  is  given  by 


V  xV=  V"  (Ui  +  I  «  vi') 
=  (U'  +  =  Qj 


The  accent,  whenever  used,  denotes  differentiation  with  respect  to  z.  • 
Henceforth  the  term  "vorticity"  will  be  applied  only  to  the  scalar  mag¬ 
nitude  of  the  vorticity  vector  since  the  orientation  is  always  along  the 
unit -vector.  The  total  vorticity  is  then  given  by 

/ 

0  =  o'  +  v'^=U'+(|  , 

where  has  been  written  for  ^  '1^  ,  the  perturbation  vorticity. 

The  models  to  be  included  are  (a)  an  unbounded  fluid,  (b)  a  fluid 
with  a  boundary  surface  below,  and  (c)  a  fluid  with  boundary  surfaces 
below  and  above.  Hereafter  these  will  be  called  the  unbounded,  ha  1  f - 
bounded,  and  boui’^dcd  models,  respectively.  In  any  of  these,  the 
vorticity  field  -with  boundary  conditions)  provides  a  unique  determina¬ 
tion  of  the  velocity  field,  and  the  perturbation  vorticity  field  is  chosen 
as  the  principal  dependent  vari.ible. 


^  Q  (x./'.t)  ■ 

'I'he  mechanism  of  the  evolution  is  contained  in  Helmholtz'  principle 
of  individual  vorticity  conservation,  DQ/D  t  =  0  •  This  is  developed 
with  the  aid  of  I  cjs.  (1),  (<i),  and  (3)  into 

“  -  lO-uoU"*)-''-'’! 


(3) 
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re 


i“  I  I  1 

(In  treating  zonal  flow,  U  Cz)  is  replaced  by  U  C?)  '  (.Ij  whe 

F(z)  is  the  Coriolis  vorticity.  LJCZjil  then  represents  a  westerly 
current,  and  X  and  2  are  the  eastward  and  northward  coordinates, 
respectively.  ) 

We  linearize  at  this  point  by  considering  the  flow  as  zero-order  state 
and  the  perturbation  as  first  order  in  smallness. 


^  q 
at  t 


dX 


(-L)(}  +  U'>)  . 


(4) 


Consequently,  we  shall  be  able  to  consider  elements  of  an  x-depen- 
dency  spectrum  singly.  The  perturbation  vorticity  of  such  an  element 
is  given  by  the  real  part  of 


^  (z.t) 


iR-X 


(5) 


where  CH,t)  is  complex. 

The  boundary  surfaces  need  not  be  smooth  but,  for  consistency  in 
this  linear  treatment,  any  departure  from  smooth  must  neceasarily  be 
of  the  lirst  order  in  smallness.  Furthermore,  the  boundary  deforma¬ 
tions  must  be  analyzed  and  paired  element  for  element  with  the  per¬ 
turbation.  Accompanying  Eq.  (5),  the  boundaries  are  given  by  the 
real  part  of 
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(6) 


where  O  is  a  complex  amplitude  factor.  So  paired,  the  mechanism 
will  not  generate  other  wave  numbers.  Completeness  is  achieved  by 
considering  the  entire  real  positive  range  of  fe,  . 

At  a  boundary  surface  the  slope  of  the  streamline  must  be  the  same 
as  the  slope  of  the  boundary.  This  condition,  integrated  with  respect 


4 


to  X 


,  gives  the  linearized  relationship 


[1=  -IJ51 


at  the  bcundary  levels.  In  the  unbounded  and  iiaU-bounded  models,  the 
additional  requirement  that  the  perturbation  velocity-field  be  upper- 
bounded  is  imposed. 

Applied  to  the  arbitrary  element,  Eq.  (4)  becomes 


t  IJ "ci)4'C2,t)]  1  , 


(8! 


and  the  perturbation  stream-function  is  then  j^iven  by 


(9) 


together  with  the  boindary  conditions  stated  in  the  preceding  paragraph. 
Since  the  x- dependencv  factori#  out  cf  the  system,  we  shall  subsc- 

j,y 

quently  omit  the  addendum  unless  specific  reference  is  desired. 

The  btream-lunction  which  satisfies  Eq.  (9)  can  be  expressed  by  a 
Green's  function  integral* 

(zt)  =  I  Cz.t)  ,  (10) 


where 

Iczp  =  . 


(11) 


♦The  Integral  can  be  deduced  b’  quasi' phys ical  reasoning  based  on  the 
superposition  of  Rayleigh  Wav' s.  See  Holl/ 


.-JP  .P..  4r jT*..*',  p-T 


For  the  general  solution  we  add  the  solutions  of  the  reduced  equation 
and  get 


t  (2,t)  =  I  (z,t)  +  +  C_  . 


the  complex  coefficients,  2ind  C_  .  of  the  harmcnic  fields  are 

determined  by  the  boundary  condtions. 

In  the  unbounded  model  the  fields  extend  from  Z  =  -oCto  Z  =0O 
Arguing  that  the  perturbation  velocity  field  remains  finite  at  2  =  ±  OO 
leads  to 

C=C=0  ^ 


The  perturbation  stream-function  for  the  unbounded  model  is  thus  given 

by  ^ 

't'cz,t)  =  r  a,t)  ,  ( 


in  which  the  subscript  denotes  the  lower  — and  the  superscript  the  upper- 
limit  of  the  integral. 

In  the  half  -bounded  model  the  fields  extend  from  the  boundary  at 
E  =  d  to  E  =  oO  . '  The  condition  of  finiteness  leads  to  ^  =  0  • 

At  the  boundary  surface  Eq.  (7)  must  be  satisfied. 


oo  ^ 

L  -U(a)  5  , 


which  determines  U_  .  The  perturbation  stream-function  for  the 
half-bounded  model  is  thus  given  by 

+  (E.t)  -  I^(Z,t)-  U(a)5+]^Kt) 


in  the  bounded  model  the  fields  extend  from  the  boundary  at  2  =  CL  t 
the  boundary  at  2  =  b  .  At  both  boundary  levels  Eq.  (7)  must  be  satis 
fied.  This  gives  a  pair  of  simultaneous  equations  which,  when  solved 
for  and  ,  yield: 

C  =  ■  2icJ  'lEH  '  -  I  (b,t) 

+  e  -e 

-  e  -  e  L 

-Ll(as5^e"MJtb)5be-^^'  • 

The  subscripts  <X  and  b  have  been  attached  to  3  so  as  to  dif¬ 
ferentiate  between  the  complex  amplitude  of  the  corrugation  below  and 
the  complex  amplitude  of  the  corrugation  above.  The  perturbation 
stream -function  for  the  bounded  model  is  thus  given  by 

=  I  ^(1.1)  t  C.«'“  +  C_e'‘‘ 

where  C.  and  C_  are  given  by  Eqa.  (14)  and  (15),  respectively. 

Regarding  the  unbounded  model,  there  is  a  special  case  which  is  of 
interest  and  which  requires  special  treatment.  This  is  the  case  of 
periodicity  in  Z.  .  The  straight  parallel  flow  and  the  initial  distur¬ 
bance  are  periodic  in  Z.  ,  and  they  have  the  same  period  so  that  the 
periodicty  will  be  maintained.  The  perturbation  vorticity  at  all  times 
satisfies  the  relationship 


where 


H  is  the  period  and  T\.  =  ...  -2,  -1,  0,  1,  2  .  . . 

To  treat  this  problem  we  need  focus  our  attention  on  only  a  single 
layer,  from  Z  =  Q,  to  Z  =  b  ,  of  thickness  H  -  b-CL  •  We 
develop  4*  (^.t)  in  this  layer,  that  is,  where  (1=  Z  =  b  » 

follows : 


oO 


-aO 


=  I 


7\. 


b  +  <vH 


OL+riH 


air.t) 

'cIl 


Put  (^  =  S  -V  n.H  »  and,  by  virtue  of  Eq.  (17), 


=l(S.t)d5. 


Subsequently, 


'l'(z,t)  = 


b 

^  „-fel 
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-ik. 
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A  A  >  -‘-'A  A'A'.S* 
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Whereupon,  summing  the  series, 


(i.t)  =  \[  (j.t)  + 

L(b,t) 

j-g-k-H 


C.(-i.t)  -k(b-z) 

4  -bM  ‘ 

l-e 

g-k(.-a) 


These  boundary-fitting  solutions  of  the  stream-function,  Eqs.  (12), 
(13),  (16),  and  (IB),  when  introduced  in  the  tendency  equation,  Eq,  (3), 
give,  in  each  case,  a  single  all-inclusive  governing  equation.  The 
equation  is  sufficiently  complicated  to  make  the  finding  of  even  very 
special  analytic  solutions  exceedingly  difficult.  Besides,  what  is 
desired  here  is  a  solution  method  which  is  applicable  to  any  arbitrary 
initial  disturbance  superimposed  on  any  arbitrary  profile.  To  at  least 
partially  achieve  this,  we  must  use  numerical  methods;  and  we  shall 
find  that  the  expressions  we  have  developed  are  well  suited  for  numer¬ 
ical  treatment. 

3.  NUMERICAL  METHODS,  AN  INTRODUCTION 

I  ' 

Finite- region  fields  in  general  cannot  be  completely  specified  by  a 
finite  number  of  pieces  of  data.  In  fact  not  even  is  the  exact  value  at 
a  point  relayable.  However,  for  numerical  tractability,  fields  must 


be  defined  (that  is,  approximated)  by  a  finite  sat  of  rounded-off  numbers 
necessarily  accompanied  by  an  interpretation  scheme.  This  discretiza¬ 
tion  can  be  accomplished  in  characteristically  different  ways. 

The  most  common  method,  that  associated  with  direct  measurement, 
is  to  give  the  value  of  a  field,  by  a  rounded  number,  at  each  of  a  finite 
number  of  points  in  the  region  and,  in  addition,  to  specify  some  interpo¬ 
lation  scheme.  Usually  there  is  some  order  in  the  spacing  of  these 
points. 

In  another  common  method  the  field  is  expret.sedas  a  linear  combi¬ 
nation  of  a  given  finite  set  of  ar.alytic  functions.  The  value  and  deriva¬ 
tives  of  each  of  these  functions  can  be  computed,  from  its  analytic  form, 
to  any  desired  accuracy  at  any  point  in  the  region.  The  field  is  defined 
by  a  set  of  rounded  numbers  which  are  interpreted,  in  a  prescribed 
order,  as  the  combination  coefficients  of  the  analytic  fuhetions  in  the 
linear  combination. 

A  fixed  number  of  pieces  of  data  represents  only  so  much  information 
no  matter  which  method  is  used.  One  method  may  give  a  better  approxi¬ 
mation  than  another  in  particular  cases,  but  basically  one  method  is  as 
capable  as  any  other.  In  fact,  the  two  methods  are  equivalent  if  inter¬ 
polation  in  the  first  is  based  on  the  fitting  of  the  functions,  used  in  the 
second,  to  the  data  at  the  given  points.  The  choice  of  method  rests  on 
peculiarities  of  the  particular  problem.  If  functions  can  be  found  which 
are  indigenous  to  the  system  of  equations,  such  that  their  use  makes  the 
problem  more  tractable,  then  the  second  method  may  De  preferable. 

Generally  the  first  method,  with  a  uniform  spacing  of  points  (called 
a  grid)  is  used.  And  it  is  often  used  with  mixed  interpolations  (low- 
order  polynomials)  even  at  the  same  place  in  a  particular  field.  It  is 
up  to  numerical  analysis  tc  determine  if  such  inconsistency  is  permis¬ 
sible.  It  is  numerical  analysis,  the  investigation  of  the  method  for 
accuracy  and  econojny,  which  gives  confidence  in  the  method  and  re¬ 
sults  therefrom. 


10 


We  shall  use  the  grid  discretization.  The  equations  and  bovindary 
conditions  which  define  the  continuous  fields  must  be  transformed  so 
that  they  define  the  discrete  fields.  This  is  generally  done  by  evalu¬ 
ating  the  equations  at  each  of  the  grid  points.  Differentials  in  these 
equations  arc  approximated  by  corresponding  differences  so  that 
algebraic  equations,  called  finite-difference  equations,  result.  The 
error  introduced  by  these  approximations  is  called  truncation  error 
because  it  can  be  regarded  as  due  to  truncating  series  representations 
of  the  derivatives. 

P.  D,  Thompson*  indicated  the  merits  in  developing  the  finite- 
difference  equations  by  an  averaging  of  the  differential  equations  over 
finite  elemental  regions  referred  to  the  grid  points.  He  has  shown 
that  certain  undesirable  biases  in  the  finite-difference  equations  are 
eliminated  in  this  way. 

A  necessary  condition  which  must  be  satisfied  by  the  finite -differ¬ 
ence  equations  is  that  they  approach  the  equations  which  they  approx¬ 
imate  as  the  space  and  time  increments  tend  toward  zero.  However, 
this  is  not  sufficient  to  insure  that  the  discrete  fields  computed  from 
the  finite-difference  equations  will  approximate  the  behavior  of  the 
continuous  fields  defined  by  the  differential  equations;  the  usurper 
here  is  computational  instability,  an  apparent  "blowing  up"  of  the 
round-off  error. 

Certain  finite-difference  equations  are  referred  to  as  "marching" 
equations  because  of  their  use.  They  are  similar  to  recurrence 
formulae.  Values  at  successive  points  are  given  in  terms  of  those 
that  came  before.  The  discretization  of  initial  value  problems  al¬ 
ways  results  in  marching  equations. 

In  well-behaved  marching  equations  the  round-off  errors  are  ran¬ 
dom  and  largely  self-cancelling,  while  in  others  they  may  grow 
rapidly  and  soon  swamp  out  all  significant  digits.  What  happens  in 
the  latter  is  actually  a  systematic  growth  of  error,  due  to  a  more 
or  less  complicated  interplay  of  round-off  errors. 

^Private  communication. 


*.  •»  -s  .  • 


•  J»  A  ^  P  J»  Jt  ■  P  ^  S  «  0 


J"  .*  -JS  -J*  -  > 


An  elementary  example  of  a  marching  equation  which  is  compu¬ 
tationally  unstable  is  the  recurrence  formula  fox'  Bessel  functions  of 
increasing  order  at  any  fixed  value  of  X  » 


J  U)-J  Cx)  5  (a=l,2., . . .)  ■ 

X  ^  r-l  ^  •  I  ' 


In  linear  theory  the  mechanism  of  computational  instability  is 
relatively  simple.  Several  forms  are  recognized.  One  is  the  admis¬ 
sion  of  extraneous  solutions  which  amplify  more  rapidly  than  the 
desired  solution.  This  is  the  cause  of  instability  in  the  above  equation. 
The  Neumann  function  of  increasing  order  is  also  a  solution  of  this 
recurrence  formula.  The  least  round-off  error  admits  the  Neumann 
function  into  the  marching  equation  and,  as  the  Bessel  component 
decreases,  the  Neumann  component  amplifies  rapidly.  In  addition,  at 
each  step,  new  round-off  error  changes  the  composition  of  the  solution. 

This  recurrence  formula  is  thus  inadequate  to  compute  more  than 
just  a  few  Bessel  terms.  On  the  other  hand,  as  a  recurrence  formula 
for  Neumann  functione,  the  above  is  computationally  stable.  Now  the 
Bessel  functions  are  the  extraneous  solutions,  but  these  dampen  and 
are  of  less  importance  than  the  purely  random  round-off  errors.  We 
shall  again  encounter  this  problem  of  extraneous  solutions  later  on. 

Some  problems  in  mathematical  phys'.cs  become  tractable  by  only 
a  partial  discretization.  We  shall  also  encounter  an  example  of  this 
situation. 

Discretization  transforms  fields  into  vectors.  Each  piece  of 
specifying  data  can  be  interpreted  as  a  component  of  this  vector  in 
some  prescribed  order.  In  a  grid  discretization  the  components  are 
given  by  the  values  at  the  grid  points.  Hence  a  field  which  is  unknown 
becomes  a  vector  with  as  many  unknown  components  as  there  are 
points  in  the  field.  And  when  the  differential  equation  which  deter¬ 
mines  this  field  is  transformed,  we  get  one  algebraic  equation  per 
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point — thus  the  same  determinateneas  is  maintained.  The  resulting 
system  of  algebraic  equations  can  generally  be  written  in  matrix 
notation,  a  form  v/hich  is  also  suitable  for  numerical  analysis, 

4.  THE  SPACE  DISCRETIZATION 


In  this  section  we  carry  out  the  discretization  of  the  2  dependency, 
the  t  -dependency  being  left  continuous.  The  Z  continuum  is  re¬ 
placed  by  the  points  2  =  Mk.M  =  0,  d:l,  :fe2,...  where  K  is  a 
fixed  increment.  The  integer  M  will  be  called  the  "address”  of  the 
point  2  =  MK  ,  and  will  be  used  as  a  labeling  subscript.  Our 


principal  dependent  field, 

Q  (t)  being  its  M  component. 

'  M 

Applied  at  the  arbitrary  point 
Eq.  (8),  tells  us  that 


becomes  a  vector. 


0  (t) 

i 


-  Mk, 


our  governing  equation. 


±  q  (t)  =  til  j-U  q 
dt  L 


Ct)  + 


(19) 


We  have  developed  four  expressions  for  4^  (2,t)  depending  on  the 
boundary  condtions.  We  are  now  going  to  develop  the  corresponding 
expressions  for  (t)  ,  to  be  given  in  terms  of 

component.  u 


We  begin  this  transformation  by  evaluating  the  integral 
at  the  arbitrary  point  Z  =  .  The  grid  is  placed  so  that  the  level 

The  fixed 


a.* 

coincides  with 


Z  =  Mh  . 

2.  ■=  (L  coincides  with  the  point  whose  address  is  M 
increment,  K,  ,  is  chosen  so  that  the  level  2-  =  b 
another  grid  point,  whose  address  shall  be  ,  and  so  that  the 

number  of  points  in  between  captures  the  desired  amount  of  detail. 

By  expressing  the  integral  as  the  sum  of  two  parts, 


.Mb 

I  (MKi)  =  I  (Mh,t)  +  I  (Mk,t) ) 
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we  can  remove  the  absolute  value  notation  in  the  integrands. 

We  shall  show  the  development  of  one  of  these  (the  second  part) 
and  proceed  as  follows: 


MK  J  “ck 


Mb-M-1 

=  I 

^  [M  h 


Mh 

M  + 1 1]  h 


-IK 


Next  replace  the  dummy  variable  ^  by  S  where  ^  S. 

W  e  get  ^ 

.M.h  ^ 

lx 


Mw-M-1  '' 

T-'b^  ..m 

1 

Mh 


rr^srO 


ci5 


) 


(20) 


k 

where  iX  -  6  .At  this  point  we  must  specify  the  interpolation 

scheme.  For  ^(S;t)  ,  in  the  interval  0  =  S  S  K  where 

varies  from  q  (t)  to  Cj  (t)  ,  )  '*'®  introduce 

4M  +  m  TMtm+l 


qCs.t)=  q  (i)  +  -^fq  Ct)  -  q  (t)  1  +  £  (21) 

4  4’M-*rr\  h  4M+fn.J  ' 

where  £  is  the  so-called  truncation  error  in  this  expression.  The 
advantage  in  using  linear  interpolation  will  be  seen  in  the  resulting 
simplicity  of  the  operational  matrix,  particularly  when  boundaries  are 
involved. 

If  (2,t)  is  analytic  in  2  ,  we  can  develop  a  series  expression 
for  the  truncation  error.  We  may  write 


Cs)  = 


(q') 


S  -4- 


(q”)  3%! 

4  M  +  tn 


(22) 
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which  i8  a  Maclaurin  series  expansion  about  Z  =  Successive 

derivatives  can  be  eliminated  by  evaluating  Eq.  (22)  at  other  grid 
points.  We  eliminate  the  first  derivative  by  substituting  for  it  from 


and  get 


i  +  d'), 


s"- 

'h+^X! 


04) 


This  gives 


=  (f), 


S 

Ts 


h(i4)"d1 

1  ®  T  M- 


4(1- 


(23) 


Proceeding  with  evaluating  our  integral,  we  introduce  Eq.  (21) 
into  Eq.  (20)  which  we  develop  as  follows: 


o\=0 

fn.«0 
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Mb-M-1 


tY\.=  0 


where 


-h.K 


K  =  ;  K  =  , 

Kk.^  ^  ^  ^ 


(24) 


and 


1 

j 


£  e  t 


k: 

II 


H+nv 


terms  of  K'lqkef 
order  in  K 


In  numerical  analyaia,  the  , value  of  judging  truncation  error  on  the 
basis  of  a  few  terms  of  the  erroil  series  is  dubious  and  this  practice  can 
be  carried  too  far.  At  the  stage  where  truncation  error  becomes 
troublesome,  the  validity  of  such  an  approximation  becomes  questionable. 
Furthermore,  there  is  a  tendency  for  such  approximations  to  give  a 
false  sense  of  security  in  higher-order  interpolation  schemes,  whereas, 
in  practice  the  arrival  of  intol8ra.ble  gross-truncation-error  in  compu¬ 
tations  is  probably  little  delayed  by  such  schemes.  More  consideration 
should  be  given  to  the  use  of  tighter  grids  ^that  is,  smaller  K,^  which 
results  in  a  more  significant  actual  reduction  of  truncation  error. 

In  making  real  predictions  from  initial  data,  some  investigators 
object  to  the  use  of  tighter  grids.  They  feel  the  number  of  grid  points 
should  not  exceed  the  number  of  initial  observations  of  the  field.  They 
also  feel  that  computations  should  be  carried  out  with  no  more  than  the 
number  of  significant  digits  in  the  observations.  We  believe  these 
notions  are  incorrect. 

It  is  important  to  distinguish  between  the  types  of  error  with  which 
we  are  dealing.  The  initial  error  in  the  fields  is  due  to  limited  ob¬ 
servations.  Truncation  error  is  due  to  the  disparity  between  finite- 
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diffe  ence  equations  and  differential  equations.  Because  we  are  trying 
to  duplicate  a  differential  process  it  behooves  us  to  use  as  tight  a  grid 
as  is  economically  feasible;  and,  enough  digits  to  give  a  significant 
difference  between  neighbouring  grid  points  should  be  carried. 

In  the  absence  of  axiy  subjective  analysis  skill,  the  initial  values 
of  the  fields  at  the  grid  points  are  computed  from  some  objective 
interpolation  scheme.  Extra  grid  points  thus  do  not  add  to  the  specifi¬ 
cation  of  the  fields  at  the  initial  moment  but  at  moments  after  the 
initial.  The  values  at  the  extra  grid  points  are  not  redundant  and  play 
a  significant  role  in  holding  down  truncation  error. 

If  it  could  be  shown  that  physical  differentia!  processes  are  trans¬ 
formable  into  algebraic  relationships  between  values  at  grid  points 
without  differentially  dependent  truncation  error,  then  the  comments 
of  the  preceding  two  paragraphs  could  be  ruled  out. 

Returning  from  this  digression  to  our  problem,  we  accept  the 
truncation  error  by  dropping  £  ,  This  error  can  be  estimated  by 
vi'iually  superimposing  linear  interpolation  on  c.  plot  of  Cj  vs.  z  . 

In  our  development  we  have  reached 


^  z  K%  t  It  Z  K"q 

In  the  second  summation,  put  TfTl  =  TA.  -  1 


Mb'M-l 


Since  m  is  only  a  dummy,  we  can  drop  the  bar,  that  is,  we  can 
replace  iTt  by  to.  .  If  we  now  express  the  first  term  of  the  first 
summation  separately  and  add  and  subtract  an  extra  term  on  its  end, 
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1>«)  ■  -It 


then  we  can  recombine  the  two  summations.  Thus 


Mh-M 


W(Mh,0  =  4  K  k  K  , 


K,  =  K,  +  4.  e  -l]  . 

K  I  J/ 


By  a  similar  treatment  we  find 


I  (M- 1)  :.-  Jh  Q  4  !ii  ) '  K”'  "‘q 

tM  +  -ZK  £-  -}M-m  -in  iMa 

Adding  the  two  parts  of  the  integral  completes  its  evaluation  at  the 
arbitrary  point  2  =  MK.  . 


MfM  ^ 


M  +  m 


+  -K,K%”  -K.Cq  ] 

rnO  '  A  '  Aj 
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The  two  other  integrals  which  need  be  evaluated  can  readily  be 
found  from  this  one.  The  first  .'esults  wh^n  Mb  goes  to  infinity. 
Eq.  (26)  becomes 


CO 


The  aecond  results  when  in  addition  we  let  go  lo  minus  infinity. 


.  oo 


oo 


I  =  :jr 

cR, 


•  oo 


,m 


+  K3Z 


m 


TtV  »  1 


•m 


We  may  now  readily  complete  the  four  expressions  for 
corresponding  to  the  different  boundary  conditions.  We  already  h.^ve 
it  for  the  unbounded  model.  According  to  Eq.  (12)  page  6,  it  is 
given  by  Eq.(28): 


VIVaL 


Th 


00 


+  K3  Z  To 

m-l  tM- 


m 


At  this  point  we  should  make  some  remark  as  to  what  is  going  to 
be  done  with  the  infinite  series.  In  computing,  these  must  of  course 
be  terminated;  and,  in  practice,  one  can  at  most  compute  the  evolu¬ 
tion  throughout  a  finite  layer  of  the  fluid  (with  the  exception  of  the 


z-periodic  case).  This  necessarily  restricts  the  type  of  profiles  that 
can  be  handled  with  this  method  if  there  is  to  be  any  standard  of  accu¬ 
racy.  Since  the  stream  function  at  any  l«vel  in  the  layer  is  to  be  com¬ 
puted  only  from  the  vorticity  in  the  layer,  it  is  required  that  the 
ignored  vorticity,  ,  outside  of  the  lay  we  are  considering  is  —  and 
is  expected  to  remain — negligible  in  comparison  with  the  vorticity  in 
the  layer.  The  method  is  valid  if  this  is  initially  so,  and  will  remain 
valid  the  longer  the  smaller  U  is  outside  the  layer.  Fortunately 
the  condition  can  oe  checked  quite  simply.  The  vorticity  at  the  end 
points  of  the  layer  may  be  regarded  as  an  indicator  of  the  magnitude  of 
the  neglected  outside  vorticity.  If  these  end  vorticities  become  ap¬ 
preciable,  again  by  comparison,  then  the  method  begins  to  break  down 
as  increasing  error  flows  in  from  the  outside. 

In  the  absence  c>f  boundaries,  tht  addresses  of  the  extreme  points 
of  the  layer  sliall  be  given  by  P  for  the  upper  and  for  tne 

lower.  The  upper  limits  of  the  summations  are  thus  replaced  by 


M  - 

'  'u 

M 

f-y 

for  oO  and  by  1 

1 

r 

c 

-  OO 

.  We  shall  also  at  this 

point 

introduce  the  notation 

M,  •  M' 

M 

(“IhI  - 

7  K  q 

tn^  1 

z_ 

f\  q  (q  ) 

tM-m  \ 

(30) 

where  the  subscripts  u.  and  1.  indicate  that  the  weighted  summation 
is  to  be  taken  over  all  the  upper  grid  points  and  over  all  the  lower  grid 
points,  respectively. 

The  stream  function  of  the  unbounded  model  is  given,  with  this 
notation,  by 

f  ,  I 

Q  's  are  functions  of  time, 

1m 


in  which  it  i  understood  that  the 


We  now  move  on  to  the  stream  function  of  the  half-bounded  model 
as  given  by  Eq.  (13)  page  6  .  The  discretization  is  accomplished  with 
the  help  of  Eq.  (27)  for  the  integrals.  Substituting  as  before  Z.  =  |V1K  , 
a  =  and  labeling  the  uppermost  point  by  ,  we  find  that 


wherein  the  term  l\  appears  twice  and  cancels.  Since  the 

summations  extend  over  all  the  grid  points  above  or  ail  the  grid 
points  below,  we  again  adopt  the  notation,  Eq,  (30).  Hence,  for  the 
half-bounded  model. 


+ 


~  -Z.il U{a)5  +  Z-K.a  + 

[  I 


In  the  same  way,  with  the  help  of  Eq,  (26),  the  stream  function 
of  the  bounded  model,  as  given  by  Eq.  (16),  page  7,  and  the  stream 
function  of  the  z-periodic  case,  as  given  by  Eq.  (18),  are  discretized. 


If  we  introduce  generalized  star  functions,  and  ,  all  four 
expressions  for  the  stream  function  can  be  given  by  one  expression: 

t,‘"  -  ^  4„  *  K3  +  1<3  (  33) 

3  h  .  3  *.  K“'"‘  . 

Lk.  ^  Zk.  c 

The  corresponding  star  functions  are  given  as  follows: 

For  the  unbounded  model, 

=  Q  . 

For  rhe  half-bounded  model, 

*1=0 

I'or  the  bounded  model, 


(34) 

(35) 
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The  explanation  of  the  new  notations  follows: 


akk)  +  h''' - 1) 


Kt  =  r  kV (K""-  K"1 

5^=  -^tiUd)  o/Kj 

A  =  K,  r  [  Uw  ?  k  r*  5,  -  U(b)  Zk  K"‘  5,  ]  /Kj 
4 = -  K,k"'  [  Ufei)  Zk  K"‘  S,-  U(b)  Zk  r  SbJ/Kj 

Substitution  of  Eq.  (33)  into  Eq.  (19)  yields  the  tendency  of  the 
arbitrary  component  ' 

^  %  ■■=  -“-[-“mHm  ''M{(q«>a+  (Oij 

-V^¥r^h  -  \  5 


where 


/*„=  ;  v„  =  U:K3A 


It  should  be  apparent  that  these  formulae  have  been  designed  for 
simplicity  in  computing.  For  a  given  problem,  the  S's  and  K's  are 
all  constants;  and  and  are  components  of  constant  vectors. 

The  5  and  the  star  functions  are  functions  of  time.  The  star 

functions  are  computed  at  any  instant  from  the  constants  and  from 

^  are  independent  of  M  . 

There  is  an  all-important  feature  on  which  rests  the  simplicity 
and  economy  of  these  formulae  for  computational  purposes.  As  can  be 
seen  from  Eq.  (30); 


^  Im-i  ■ 


(44) 


It  is  these  recurrence  formulae  which  make  it  possible  to  compute 
the  tendency,  at  all  grid  points,  by  making  at  most  only  five  computing 
traverses  of  the  grid  points,  no  matter  how  many  points  there  are  in  all.  The 
procedure  is  as  follows: 

a.  Beginning  with  the  top  point  and  progressing  downward,  the  con¬ 
tribution  of  the  term  containing  (  )i*.  computed  at  each  point 

with  the  help  of  Eq,  (44). 

b.  Beginning  at  the  bottom,  the  contribution  of  the  term  containing 
(c^^)^  is  computed  at  each  point  also  making  use  of  Eq.  (44).  Upon 
completion  of  these  first  two  traverses  the  two  star  functions  are 
evaluated. 

c.  Then,  beginning  at  the  top,  the  contribution  of  the  term  containing 

Kf'^i  -  M 

is  computed  at  each  point. 

d.  Beginning  at  the  bottom,  the  contribution  of  the  term  containing 

computed  at  each  point.  These  latter  two  traverses 
require  only  a  multiplication  by  K  in  progressing  with  the  factor  from 
one  point  to  the  next. 


24 


O  'J  V  ■  /  f .  W“.  r,  e-.  •“«  r# 


w 


) 


I 


) 

I 


5 


1.. 


e.  The  contribution  of  the  term  containing  may  be  obtained  by  a 

separate  traverse,  or  may  be  combined  with  one  of  the  other  traverses. 

At  this  stage,  we  may  ask  about  the  applicability  of  this  investiga- 

.11 

tion  method.  For  one  thing,  U  must  be  continuous.  Should  we 
desire  to  investigate  particular  profiles  in  which  this  is  not  the  case, 
then  we  must  either  use  a  modified  procedure  or  we  can  make  an 
approximation  to  the  flow  profile  by  fitting  to  it  a  curve  which  has  a 
continuous  second  derivative. 

The  perturbatsjn  vorticity  distribution  in  c  is  being  approximated 
by  a  linear  interpolation  between  grid  points.  The  magnitude  of  the 
apace  increment,  K  .  must  be  sufficiently  small  so  that  linear  inter¬ 
polation  will  closely  approximate  the  initial  distribution.  Whether 
the  interpolation  will  yield  a  good  approximation  at  later  times  de¬ 
pends  also  on  the  smoothness  of  UlZ)  and  U  (2^.  Hence,  h  must 
be  chosen  sufficiently  small  to  also  capture  the  detail  of  these  functions. 
Even  with  these  precautions  the  linear  interpolation  of  |  during  its 
evolution  may  become  completely  inadequate.  This  is  called  the  arrival 
of  the  intolerable  gross-trimcation- error. 

5.  EIGENSOLUTIONS  AND  COMPUTATIONAL  STABILITY 

For  purposes  of  numerical  analysis  the  governing  set  of  linear 
equations  which  is  represented  by  £q.  (42)  is  written  in  matrix  notation; 


q(t)  =  -  L 
dt  1 


C  q(0  d 


(45) 


The  coefficient  matrix,  C  ,  has  only  r«il  constant  elements.  The  con- 

I 

stant  vector,  q  ,  arises  from  the  deformation  of  the  boundary  sur¬ 
faces  and  may  be  complex.  The  presence  of  the  boundary  surfaces,  but 

not  their  deformation,  affects  the  elements  of  C  . 

df 

The  evolution  of  initial  perturbations  can  be  determined  without 
discretizing  the  t  continuum.  Partial  discretization  has  made  the  prob¬ 
lem  tractable  but  the  work  is  prohibitive  unless  has  few 
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components.  The  solution  method  is  based  on  the  expansion  of  the 
initial  field  into  a  linear  combination  of  eigensolutions. 

A  particular  solution  of  Eq.  (45)  is  the  steady  state 

-1 

=  -  C  d  .  (46) 

This  may  not  be  the  most  general  steady  state  with  wave  number  L  . 

For  the  complete  solution  of  Eq.  (45)  we  must  add  to  the 

general  solution  of  the  reduced  equation, 

i  q(t)  -  -  L  C  q(A)  ■  (47) 

dt  i  ~  i 

Introducing  the  eigensolution 

C^(t)  =  ,  (48) 

into  Eq.  (47)  leads  to 

"  Xq  =  ~  Cq 

T  4 

This  is  equivalent  to 

( 9  ■  V  9  ’ 

where  is  the  identity  (or  unit)  matrix. 

This  system  of  linear  eqiiations,  Eq.  (49).  being  homogeneous 
is  overprescribed  (thus  cannot  have  any  solution)  unless  the  deter¬ 
minant  of  the  coefficient  matrix  vanishes: 

I  C  -  'Xl  I  -  0  •  (50) 

This  is  a  polynomial  in  X.  and  is  called  the  characteristic  polynomial 
of  the  matrix  C  .  Its  order  is  that  of  the  matrix.  Thus  it  has  in 
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general  as  many  different  roots  as  we  have  grid  points.  These  are 
called  the  latent  roots  of  C  and  are  the  eigenvalues  which  we  seek. 
Corresponding  to  each  latent  root  there  is  a  non-zero  vector 

which  satisfies 


(51) 


These  are  the  latent  vectors  of  L  and  are  the  eigenvectors  which 
we  seek.  They  are  indeterminate  to  the  extent  of  an  arbitrary  factor. 
Thus  we  may  add  that  they  be  normalized,  that  is,  |  |  =  1  . 

If  Eq.  (50)  admits  X  =  0  as  a  solution,  then  the  corresponding 
eigenvector  can  be  combined  with  q  to  form  a  more  general  steady 


state.  ^ 

Because,  in  general,  there  are  as  many  different  eigensolutions  as 
there  are  grid  points,  the  set  of  eigensolutions  is  complete.  That  is, 
by  a  proper  choice  of  complex  weighting  coefficients,  a  linear  combina¬ 
tion  of  the  eigensolutions  added  to  can  be  made  to  fit  any  initial 

vector  and  thus  will  give  its''subsequent  evolution. 

Because  C  has  all  real  elements,  its  latent  roots  are  either 
real  or  occur  in  complex  conjugate  pairs.  If,  as  a  special  case,  C 
is  symmetric  aa  well,  then  all  its  roots  must  be  real.  However,  the 
roots  of  an  asymmetric  matrix  may  also  all  be  real. 

For  a  real  latent  root  it  is  clear  that  the  corresponding  latent 
vector  can  always  be  chosen  so  as  to  be  real.  If  not  so  chosen,  its  real 
and  imaginary  parts  must  be  parallel. 

If  ^  corresponds  to  one  latent  root  of  a  complex 

conjugate  pair  then  %  -  A  Q  —  LBq  corresponds  to  the  other, 

/\  and  Bo  here  being  real  and  not  parallel.  These  statements 
can  be  verified  by  substitution  in  Eq.  (51). 

The  above  remarks  have  a  bearing  on  the  "tilt"  of  the  eigensolutions. 
To  introduce  the  tilt  concept  we  must  first  return  our  X  dependency. 
The  eigensolution  is  then  given  by 
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<}  (x.t)  =  Ri'.  e 

The  substitution  of  Q  =  Aq  +  iBo  and  X  =  LI  +  IV.  where 

A  Q  I  ~  ^  ' 

Aq'  IJq  »  p-  V  are  real,  leads  to 

c\  [  A^  cotL  (kx  -  jit)  '  Bo  ^  (1^^  'M)  ] 


(52) 


=  Po  (b->At  -  0o) 


The  components  of  Uq  ,  the  amplitude  vector,  and 
phase-angle  vector,  are  given  by 


D. 

<Pc 


(a:. 


,  the 


If  varies  with  M  (t  liat  is,  with  the  height  Z  )  we  say  that 

the  wave  "tilts,  "  The  shape  of  the  nodal  "line"  is  given  by  the  nodal 
vector  X  =  /»L  • 

The  eigensolution  which  has  a  real  eigenvalue  (a  real  latent  root) 
yields  a  neutral  wave  which  has  no  tilt.  The  eigensolution  which  has  a 
complex  eigenvalue  (one  of  a  complex  conjugate  pair  of  latent  roots) 
yields  a  tilting  wave. 

A  complex  conjugate  pair  of  eigenvalues  yields  one  amplifying 
wave  and  one  damping  wave.  These  have  parallel  amplitude  vectors  and 
equal  but  opposite  tilt.  Where  the  nodal  line  of  one  tilts  forward,  the 
nodal  line  of  the  other  tilts  backward  by  the  same  amount. 

The  dynamic  stability  properties  of  the  flow  profile  a  e  revealed 
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by  the  eigensolutions.  In  actual  fluids,  flow  of  smaller  scale  (down  to 
and  including  Brownian  motion)  will  generate  them.  Thus,  if  Eq.  (50) 
admits  any  complex  conjugate  pair  of  latent  roots,  the  perturbation 
will  grow  indefinitely  as  it  contains  an  exponentially  amplifying  eigen- 
solution.  The  profile  must  then  be  considered  a  transient  state  as  it 
is  dynamically  unstable. 

The  eigensolution  synthesis  method  of  solving  an  initial-value 
problem  presents  a  formidable  task  which  might  be  entertained  if  we 
are  dealing  with  a  small  number  of  components.  It  is  more  feasible 
to  complete  the  discretization.  Discretizing  the  time  dependency 
results  in  a  marching  equation. 

In  proceeding  to  discretize,  one  discovers  that  for  our  system 
there  are  a  number  of  finite-difference  analogues  which  satisfy  the 


necessary  condition  for  validity  mentioned  on  page  11.  The  choice 

among  these  is  made  according  to  which  of  the  analogues  preserves 

the  nature  of  the  time-continuous  system  as  revealed  by  its  eigen- 

solutions.  We  shall  examine  some  of  the  arialogues  in  this  manner 

4 

which  has  been  expounded  by  Hyman  and  others. 

The  t  continuum  is  replaced  by  the  discrete  values  t  -  NT  where 
N  is  an  integer  and  X  is  the  time  increment.  The  integer  N 


will  be  used  as  a  labeling  subscript.  The  continuously  varying  vector 
(t)  is  thus  replaced  by  discrete  values;  at  the  time  t 

^  The  resulting  marching  equation  also  has  characteristic  solutions. 


We  shall  call  these  the  X  eigensolutions.  They  have  the  form 


(53) 


whereas  the  eigensolutions  evaluated  at  t  =  NT  are  given  by 


^  CnT)  =  (e 


(54) 
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We  shall  see  that  the  eigenvectors  are  unchanged  by  the  discretization. 
Hence  the  change  suffered  by  a  particular  eigensolution  is  entirely  re¬ 
vealed  by  the  disparity  between  e  and  Y  . 

All  the  valid  analogues  of  Eq.  (45),  page  25/  liave  the  same  partic¬ 
ular  solution  as  has  Eq.  (45)  as  given  by  Eq.  (46).  For  this  reason  we 
shall  refer  the  analogues  directly  to  the  reduced  system,  Eq.  (47), 
thus  avoiding  the  repetitious  reduction  in  each  case. 

The  first  analogue  we  shall  examine  results  when  a  forward  first- 
order  difference  is  introduced  for  the  derivative.  The  system, 

Eq.  (47),  becomes 


Substitution  of  Eq.  (53)  results  in 


Hence 


zJf 

X 


(55) 


which  restricts  V  .  Comparison  with  Eq.  (50)  reveals  that 
( 1  -  ry  (.T  can  be  identified  with  X  .  That  is,  corre¬ 
sponding  to  each  eigensolution  with  eigenvalue  satisfying 

Eq.  (50)  we  have  a  TT  eigensolution  whose  't  eigenvalue,  , 
is  given  by 

. 

Also  revealed  by  Eq.  (55)  is  that  the  "X  eigenvectors  are  the  same 
as  the  corresponding  eigenvectors. 
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Corresponding  to  a  real  X  .  Iri  =  1  ^  (txr  ,  that  is,  the 

magnitude  of  IT  is  greater  than  one.  This  means  that  all  neutral 
waves  are  converted  into  amplifying  waves  by  the  discretization.  We 
need  investigate  no  further;  this  analogue  is  computationally  linstable. 

We  next  examine  the  centered  second-order  analogue  which  has 
received  such  wide  use  for  a  first -order  system  that  it  can  be  called 
the  "conventional"  analogue.  For  Eq.  (47)  it  takes  the  form 


0.  -  q 

^  j'M  +  l  +N-1 


Substitution  of  Eq.  (5  3)  results  in 


1 

_ 

iirr.  4  ) 

=  0  . 


Thus,  according  to  Eq.  (50)  the  X  's  are  given  by 


(56) 


that  is 


r  , -I  'A 

±  [l-(tXa)]  -  '■■rKn 

Each  X  gives  rise  to  two  values  of  for  the  same  eigen¬ 
vector - twice  as  many  'X  eigensolutions  as  there  are  eigensolu- 

tions!  This  multiplicity  is  needed  for  completeness  because  both 
the  initial-value  vector  and  the  vector  at  K)  =  1  must  be  accom¬ 
modated.  Both  these  vectors  must  be  specified  before  one  can  begin 
to  march  with  Eq,  (56). 

This  analogue  hardly  preserves  the  nature  of  the  first-order 
tendency  system  which  it  purports  to  approximate.  It  actually  has  a 
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nature  corresponding  to  a  second-order  tendency  equation. 

Comparison  with  6  ''  shows  that  the  roots  given  by  the  plus 

sign. 

r  2l 

=  +  ['  -  C-X,,)_  (57) 

are  the  proper  roots.  But  we  must  attach  a  condition  for  the  com¬ 
putational  stability  of  this  set  of  roots.  For  a  real  (neutral 

wave) . 

I  -  t  -  (tX/  +  (tx/  -  I 

provided  the  square  root  of  Eq.  (57)  is  real,  that  is,  provided 

t  ^  l/|x|  (58) 

According  to  Rayleigh's  theorem,  neutral  waves  move  with  the 
speed  of  the  basic  current  at  some  level  in  the  flow.  *  The  iransla- 

y 

tion  speed  of  the  neutral  wave  is  given  by  Xy^k  »  hence 

|x|  i  tijul  _  . 

mojf.  '  n\OJ/ 

Introduced  into  Eq.  (58),  this  yields  the  possibly  stronger  condition 
(that  is,  sufficient  but  perhaps  not  necessary): 

-tsi/xluL^ 

For  a  given  T  this  defines  a  critical  wave  number; 

k  =  I A  I  u  . 

c  /  I  may 
♦See,  e.  g.  ,  Garcia^,  p.  90, 
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Neutral  waves  with  a  longer  wave  length  (  )  are  preserved 

neutral  but  those  with  a  shorter  wave  length  (  It,  >  I 
computationally  amplified. 

The  presence  of  the  extraneous  roots  is  far  more  cause  for 
alarm.  Were  it  not  for  round-off,  ^  could  be  chosen  in  complete 
accord  with  and  the  proper  %  eigensolutions.  Unfortunately, 

the  extraneous  T  eigensolutions  do  enter  into  the  computing  as  a 
result  of  round-off,  or  more  rapidly  by  an  incompatible  first-time- 
step  (for  example,  a  forward  step). 

Corresponding  ro  a  propagating  neutral  wave  (  real)  the 
extraneous  T  eigensolution  introduces  an  additional  neutral  -vave, 

I  j  =r  I  ,  travelling  in  the  opposite  direction. 

Corresponding  to  an  amplifying  wave  (one  of  a  complex  conjugate 
pair)  the  extraneous  wave  has  the  same  tilt  but  dampens.  Corres¬ 
ponding  to  a  damping  wave  (the  ether  of  this  conjugate  pair)  the 
extraneous  wave  lias  the  same  tilt  but  amplifies. 

It  is  fortuitous  that  the  latent  roots  in  our  case  occur  in  complex 
conjugate  pairs;  otherwise  the  "conventional”  analogue  could  introduce 
amplifying  waves  (from  damping  waves)  v'hcre  none  witn  tluit  wave 
length  should  be  present.  This  could  happen  in  other  problems.  As 
it  is  here,  the  main  nuisance  value  of  the  extraneous  'X  eigen- 
aolutions  is  the  distortion  they  can  cause  in  the  amplifying  configura¬ 
tions. 

The  fact  that  the  extraneous  unstable  component  is  oscillating 
{  negative),  that  is,  changes  sign  with  each  step,  may  make  this 
trouble  readily  detectable.  The  conventional  analogue  is  further 
criticized  in  Section  10. 

Extraneous  X  eigensolutions  appear  whenever  a  finite-differ¬ 
ence  analogue  relates  more  instances  in  time  than  are  permitted  by 
the  order  of  the  differential  equation.  A  first-order  differential 
equation  by  its  nature  permits  only  a  two-point  analogue. 


A  number  of  analogues  which  are  composed  of  two  operations  have 
been  prepared.  These  readily  yield  to  analysis.  One  particular  group 
has  in  common  the  second  operation 

where  estimate  of  the  value  being  sought  and  wlxich 

may  be  biTsed  on  one  of  several  different  first  operations. 

The  estimate  may  come  from  a  simple  linear  extrapolation  of 
the  vector  itself: 


q  ^  C|  ~  q 

rN  +M'I 


(60) 


Substitution  of  Eq.  (60)  in  Eq.  (59)  reveals  the  true  form  of  this 
analogue: 


V  vjA  • 


It  is  a  three-point  formula  which  will  have  extraneous  X  eigen- 
solutions  Substitution  of  Eq.  (53)  leads  to 


0 


Hence,  according  to  Eq.  (50),  page  26,  the  corresponding  Y'b  are 
given  by 


(61) 
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I 

r: 

I , 


A 


-i 


Another  possible  combination  results  when  the  estimate  comes 
from  a  forward  time-step: 


Combined  with  £q.  (59)>  the  result  is 


[%.r  iM]A  =  - 


Introduction  of  Eq.  (53)  leads  to 


itC  ~  ^  =  5  > 


and  the  condition  on  V  is 

/ 


[t‘C‘+-itc]-  (i-r)I 


=  0  . 


(62) 
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According  to  a  useful  theorem*  of  matrix  algebra,  the  latent 
roots  of  the  matrix  | 
and  the  latent  vectors 

^  “  ^TV  •  (64) 

A  third  combination  results  when  the  estimate  for  £q.  (59) 
comes  from  the  conventional  analogue: 


'X'^  +  i-T’C  j  are  given  by 

are  those  of  C  .  Hence, 


.  7 

*See,  e.  g.  ,  Milne,  ,  p.  163. 


.S  .\\\’  -w 
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This  leads  to 


None  of  these  schemes  is  satisfactory.  The  rea4er  may  complete 
the  critique  by  mapping  their  respective  transforms,  Eqa.  (61),  (64), 
and  (65), 

We  now  move  on  to  the  most  natural  finite-difference  analogue  of 
our  first-order  system.  This  is  the  “centered"  first-order  scheme: 


The  significant  feature  is  that  the  difference  is  related  to  the  tendency 
evaluated  for  a  centered  mean-value  of  the  argument. 

Introduction  of  Eq.  (53)  leads  to 

(^“  j^l)  %=  2  • 

Hence, 

2.  I  ~  Vt\.  _  N 

li-Trv  “  ’ 

that  is, 

T  —  ^  ~  ^n.  / 

one  value  of  corresponds  to  each  same  eigenvector. 

To  simplify  the  investigation  we  introduce  XT/z.  =  +  jLoJ  , 

where  1]  and  CJ  are  real.  'iT  becomes 


This  is  to  be  compared  with 

-1\T  Z(ji  IC-2t7^ 

e  ^  e  e 


Exammation  reveals  excellent  agreement  and  no  computational 
instability.  Neutral  eigensolutions,  CC  -  0  ,  become  neutral 
'X  eigensolutions.  Damping  eigensolutions,  Co  <  0  ,  becomes 
damping  X  eigensolutions.  Amplifying  eigensolutions,  CO  >  0  i 
beoome  amplifying  X  eigensolutions. 

For  I  ^  ^  ^  quantitative  agreements  of  the 

corresponding  phase-speeds  and  the  corresponding  growth-rates 
are  very  good  indeed.  Therein  lies  the  advantage  of  using  small 
values  of  T  ,  the  time  Increment. 

Equation  (68)  can  be  used  to  compute  7^  and  CO  from  ^  for 
a  particular  eigensolution. 

This  investigation  reveals  the  centered  first-order  finite-differ¬ 
ence  analogue  to  be  the  most  natural  analogue  of  our  first-order 
system.  It  will  correctly  detect  the  dynamic  instability  of  flow  pro¬ 
files. 

This  is  the  analogue  we  shall  use. 

It  may  be  feasible  in  special  cases  to  transform  a  first-order 
system  into  a  legitimate  second-order  system.  This  is  necessarily 
accompanied  by  a  corresponding  reduction  in  the  number  of  variables. 

A  second-order  derivative  has  an  easy-to-use  analogue.  Further 
more,  the  conventional  analogue  may  be  used  for  a  first-order  deriva 
tive  when  it  appears  in  a  second-order  equation. 


The  desired  transformation  is  not  accomplished  by  merely  raising 
the  order  of  the  system  without  losing  variables  in-the  process.  If  the 
same  number  of  variables  are  maintained  in  the  resulting  higher-order 
system,  then  the  nature  of  the  system  has  thereby  been  corrupted  and 
extraneous  eigensolutions  have  entered.  The  extraneous  eigensolutions 
are  carried  over  into  T  eigensolutions  and  they  would  subsequently 
develop  in  the  computations. 

The  variables  lost  in  the  desired  transformation  manifest  themselves 
as  initial  tendencies  and  are  related  to  the  tendencies  during  the  evolution. 

Such  a  transformation  is  particularly  simple  in  our  system  because 
C  is  real.  The  differentiation  of  Eq.  (45),  page  25,  and  the  subse- 
quent  substitution  of  Eq.  (45)  for  the  first  derivative  which  appears  on 
the  right,  leads  to: 

=-C  [C<|tt)  +  d]  .  (6„ 

The  real  part  of  Eq.  (69)  and  also  the  imaginary  part  thereof  form 
two  complete  systems,  each  having  half  the  number  of  variables  of  the 
foimer  system,  Eq.  (45). 

If  we  substitute  (t)  —  *  where  and  BC’t) 

are  real,  the  real  part  of  Eq.  (69)  is 

C  KCt)  dj  ,  (70) 

for  which  we  have  an  initial  value.  A©  ,  according  to  Eq.  (45), 
an  initial  tendency 

( ^  ACo\  =  B  +  ilju  d 

\dt~  y,  — 


4  Aft)  =  -  c 
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The  lost  variables,  the  components  of  B«)  ,  are  related  to 


at  all  times  by 


i  Act)  =  Bet)  +  itn/;  d  , 
dt  ~ 


this  relationship  being  the  real  part  of  £q.  (45). 

In  the  same  way  we  could  deal  with  BC't)  a-®  a-  principal  dependent 
variable. 

Equation  (70)  yields  the  same  particular  solution  for  as  did 

Eq.  (45),  namely,  Ac,  =  ~C  Ri.;  d  .  To  determine  the  general 
solution  we  require  a  complete  set  of  eigensolutions  of  the  reduced 
equation 

-Ft 

If  Q  €!***■  ‘  is  a  complex  solution  of  Eq.  (71)  then  the  real  and  the 
imaginary  parts  thereof  are  separately  real  eigensolutions  of  Eq.  (71). 
Substitutions  of  sa*''  in  Eq.  (71)  leads  to 


c -r  r  =0  . 

I  * 

I  ^ 

Thus,  according  to  Eq.  (50),  page  26,  ±  are  the  latent  roots 

which  correspond  to  the  latent  vector,  cf  C  .  The  roots  —  Xrv. 

do  not  add  anything  alien  or  extraneous,  they  naturally  belong.  Their 
presence  in  Eq.  (71)  verifies  that  the  complex  latent  roots  of  the  real 
matrix  Q  must  occur  in  conjugate  pairs,  or  else  the  nature  of  the 
problem  would  indeed  have  changed. 

The  finite -difference  analogue  of  Eq,  (70)  after  reduction  is 
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Tt4 

^  is  a  complex  solution  of  Eq.  (72)  then  the  real  and  the 

imaginary  parts  thereof  are  separately  real  X  eigensolutions  of  Eq. 

Yn 

^  in  Eq.  (72)  leads  to 


zr-i-r^T 
T"r  i 


=  o. 


Thus,  according  to  Eq.  (50) 


i 

TV 


which,  solved  for  7^  .  yields 


=  1  -  1 Xt,  ±  t  X  T+X*  ) 


i/i 


corresponding  to  the  eigenvector  .  The  comparison  of  these 

two  values  of  7^  with  shows  the  necessary  quali> 

tative  as  v/eli  as  excellent  quantitative  agreement.  The  analogue  is 
computatior.ally  stable. 

The  application  of  this  analogue  is  quite  straightforward.  Each 
time-step  requires  the  procedure  described  on  page  24  to  be  per¬ 
formed  twice  with  slight  modification. 

Although  it  is  highly  recommended,  we  shall  not  use  this  analogue 
of  the  cultivated  second-order  system  because  we  made  use  of  special 
features  of  our  first-order  system.  We  are  interested  in  dealing 
directly  with  first-order  systems  for  which  we  use  the  natural  cen¬ 
tered  first-order  analogue  which  we  have  analysed. 


(73) 
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6.  THE  INTEGRATING  PROCEDURE  (LINEAR  PROBLEM) 


The  centered  first-order  analogue,  • 


1 

_ 

1 

II 

1 

'Ll 

C 

-W  -J 

\/z.  +  d  - 

can  be  made  explicit  as  follows: 

/  \ 


f  ItC  q  =  2.1-  i.tC  q  -  i2Xd 


/  Z' 


(74) 


ll  +  ltC 


Zl-LtC  i-lll 


f  itC]  L 


w  ,  iZ.f  d 


In  some  cases  it  might  be  desirable  to  use  this  explicit  formula. 
However,  several  objections  may  be  raised;  (a)  the  inversion  and 
I  multiplication  involve  considerable  calculation  with  accompanying 

inaccuracies:  (b)  the  magnitude  of  the  time  step  is  frozen;  and  {c)  the 
.  resulting  coefficient  matrix  on  the  right  hand  side  is  generally  far 

more  complicated  than  the  matrix  C  <  Simplifying  features  such 
as  exhibited  by  Eq,*(44),  page  24,  are  lost. 

Because  of  Eq.  (44),  it  may  be  more  convenient  to  use  an  itera¬ 
tive  method  to  achieve  Eq.  (74).  This  method  is  particularly  inter¬ 
esting  because  of  its  applicability  to  nonlinear  systems  (see  Section  11). 

Equation  (74)  may  be  written 


(75) 
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V.  V.  f 


The  iterative  scheme  for  arriving  at  the  solution  of  this  equation  is 
revealed  by 


w 


(76) 


where  the  superscripts  in  parenthesis  label  the  successive  approxi¬ 
mations  of  • 

To  analyze  the  convergence  of  this  scheme  we  set 


q  =  6  +  •,  A.  =  0  I  2. . .  .  (77) 

tn+I  XNtl  ~ 


and  thereby  define  as  the  error  of  the  fL-th  estimate.  The 

subtraction  of  Eq.  (7|)  from  Eq.  (76)  reveals 


(78) 


Hence  by  induction 


(79) 


This  scheme  is  convergent  if  0  cw  h, — >90  .  The  neces¬ 

sary  and  sufficient  condition  for  convergence  is  that  all  the  latent  roots 
of  the  matrix  ^  J  less  than  "one"  in  magnitude. 

According  to  Eq.  (50),  page  26,  the  latent  roots  of  the  matrix, 

("  =)^  given  by  ^  where  tv  =  1,  2  ....tuq 

and  is  the  order  of  the  matrix.  Its  latent  vectors  are  those  of  C  : 
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Let  us  impose  that  the  latent  roots  be  numbered  in  descending 
order  of  ma gnitude ;  and  let  us  express  .  6^®^  as  a  linear  combination 
of  the  latent  vectors: 

rv  T\.  +TL 

Thus,  according  to  Eqs.  (79)  and  (80), 


e<*>  =  Z  fe  (-i-iXXo  . 

~  n  H  2.  >71, 


This  makes  apparent  the  stated  necessary  and  sufficient  condition  for 
convergence.  Imposing  this  epndition  on  Eq.  (82)  yields 


T  <  z/ix  i 


This  is  a  definite  limitation  on  the  size  of  the  time  increment  which 
may  be  used  with  this  iterative  scheme. 

On  page  32  we  stated  that  according  to  Rayleigh 

if  X|  i^aal.  If  we  assume  tha.t  this  also  includes  complex, 

we  are  led  from  Eq.  (83)  to  the  sufficient  condition 


X  < 


z/k-l 


The  first  guess,  chosen  objectively.  Let  us 

consider  ;  nd  compare  choosing  extrapolation  (which  we  shall  indicate 
by  E:)  with  choosing  persistence  (which  we  shall  indicate  by  P:). 


P  -  =  q 

XMtl 


•N  V> -.N  V-  SI*.-  S  S  S' ■s-  S'".  S'  S-  S-  s-  s'  S'  S'  S'  s"  s'  ^  - 


This  yields 


E  •  + 1  ^  *1  “  ^ 

tH-1  .  itJ  +  l  ■ 

P :  e  =  n  -  n 

~N-VL  JM  ^N-il 

We  can  expand  as  a  weighted  sum  of  X  eigensolutions,  Eq. 

(53),  page  29, 

where  the  (X  's  are  determined  by  fitting  the  initial  vector  and 

the  f  's  are  given  by  Eq.  (67).  It  follows  that  the  choices  yield 

E  -  (f.U  = -cr^C-r/ C  . 

(85) 

P--  • 


It  is  of  interest  to  determine  which  of  the  two  cnoices  yields 
smaller  values  for  the  's.  The  ratio  of  the  members  of  Eq.  (85) 

yields 


This  magnitude  is  ^  1  for  Re;  ^  i/2,  .  Thus  Re;  1/2 

is  critical.  According  to  Eq.  (68), 


SU:T  = 


(I  ^  Cxi)  (l-Oj)-W^ 


•> 
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where 


XTy^Z.  —  77  +  L  6l)  ;  and  the  critical  value  lies  on  the  circle 


rf-  ^  (cl)  -  1/3)^  =  (.l/sf 

For  values  of  X^T/2.  falling  inside  the  circle,  extrapolation 
yields  a  smaller  value  of  than  does  persistence.  That  all 

values  of  fall  inside  the  circle  can  be  assured  by  the  suf¬ 
ficient  condition,  j  Xj  K  1/3,  that  is,  by 

X<  2/3|xJ  (86) 


In  practice,  the  incidence  of  two  successive  estimates  which  differ 
in  each  component  by  less  than  some  presct'ibed  tolerance  magnitude, 

T  ,  shall  be  accepted  as  convergence.  The  iteration  which  achieves 
this  is  denoted  by  the  subscript  C  >  that  is,  by  ,  and  the  ac- 

cepted  solution  is 

If  the  matrix  Q,  is  symmetric,  then  its  latent  roots  are  mutually 
orthogonal  and  the  analysis  of  the  error  extinction  is  rather  elementary. 
However,  this  is  too  special  a  case  to  dwell  on  and  we  shall  not  be  so 
restrictive. 

If  the  matrix  C  is  asymmetric,  then  the  magnitude  of  the  error 
vector  and  its  components  do  not  necessarily  decrease  monotonically 
until  the  slowest  decaying  component  of  the  error  latent-vector  expan¬ 
sion  dominates. 

After  some  minimum  number  of  iterations,  , 


(-  I  i 

If  ^  analysis  of  the  rate  of  convergence  can  be  based 

entirely  on  Eq.  (87)  because  what  happens  to  the  other  latent-vectors 
of  the  error  expansion  does  not  matter. 

Equation  (85)  reveals  that  in  both  cases  the  weights  of  the  error 
latent-vector  expansion  develop  with  the  "X  eigensolutions. 
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Corresponding  to  neutral  waves,  |  1  =  1  and  the  does  not  grow. 

But  in  the  case  of  dynamic  instability  |  ^  |  ^  some  I  ,  and 

the  corresponding  G\_  grows  to  progressively  increase  .  For 

this  reason  or  more  directly  because  T  may  Rave  been  chosen  in¬ 
sufficiently  small,  may  occur  before  ^  happens, 

an  effective  analysis,  with  C  asymmetric,  depends  on  all  particu¬ 
lars  of  the  problem  and  is  tantamount  to  solving  the  particular  evolution. 
We  restrict  ourselves  to  A.^  >  . 

The  reiteration  is  ended  when 


M  I 


CM  _  (A,-.) 

+N+I  tN+1 


(88) 


has  all  components  less  than 
because 


T"  in  magnitude.  This  limits  the  error 


A 


(At) 

’Ntl 


_  g(Ac-') 

~N  +  I 


(-If  x.-i) 


Up  to  this  point  we  have  not  discussed  the  effect  of  round-off  during 
the  reiteration.  Some  comments  are  now  called  for.  Round-off  causes 
spurious  variations  which  are  usually  confined  to  the  last  few  binary 
digits  of  the  components  of  during  each  iteration.  Because  of 

these  variations  the  iterations  may  not  lead  to  a  steady  value  for 
Instead  the  reiteration  may  settle  into  a  cycle  of  two  or  more  iterations. 
It  may  happen  that  during  such  a  cycle  none  of  the  successive  differ¬ 
ences  satisfies  the  tolerance,  and  convergence  will  not  occur.  This 
bothersome  phenomenon  should  not  be  mistaken  for  divergence.  It 
causes  oscillations,  not  "blow-up." 


46 


This  is  a  shortcoming  of  the  test  method  to  determine  satisiactory 
convergence.  The  probability  of  the  occurrence  of  such  oscillations  is 
reduced  by  increasing  T  ,  On  the  other  hand,  it  is  undesirable  to 
increase  the  tolerance  not  only  because  of  reduced  disparity  in  Eq.  {74) 
but  also  because  the  terms  of  the  centerv’id  first-order  analogue  will  not 

(q) 

have  been  met.  In  fact,  if  +  l  b  jing  determined  by  extrapola¬ 

tion  then  our  operation  still  relates  three  instances  in  time.  The 
tolerance  should  be  of  the  order  of  round-off  so  that  only  error  which  is 
practically  random  remains. 

Instead  of  testing  for  convergence  it  may  be  preferable  to  aimpiy 
carry  out  and  accept  the  result  of  a  fixed  number  of  iterations  for  each 
and  every  time  step.  The  number  of  iterations  may  be  determined  by 
fixing  a  maximum  error-to-wave-amplit”de  ratio. 

If  extrapolation  is  used  the  error  vector  is  given  by 


and  the  vector  itself  by 


Consequently  the  errbr-to-wave-amplitude  ratio  for  the  n-th  wave 
(n-th  latent  vector)  is 


We  set  the  upper  limit  £  on  the  magnitude  of  this  ratio  fo  all  tv  . 
If  Eq.  (86)  is  satisfied,  it  is  sufficient  that  | '^  j  ^  E  i  or 


A.  =  inv  E  Inv  }  ^ 

Both  of  these  jErv,  's  are  negative. 
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This  method  also  takes  on  more  significance  if  ^  is  symmetric, 
otherwise  it  also  has  its  shortcomings.  It  is  not  as  readily  extensible 
to  the  iterative  solution  of  nonlinear  equations  as  is  the  test  method.  In 
tne  integrations  to  follow,  the  test  method  is  used.* 

The  amount  of  computing  may  be  reduced  by  making  use  of  the  fact 
that  C  has  all  real  elements.  We  introduce  Q  =  ±  ") 

where  and  D^.  are  real,  into  Eq.  (76),  page  42.  The  real 

and  imaginary  parts  of  Eq.  (76)  then  yield 


A*!"' ,  =  A  1  [cfB  +  • 

~  M -H  ^  Ivl  I  ~  N  N+aj/  j 


(89) 


This  reveals  that  it  is  redundant  to  compute  all  indicated  successive 

estimates  for  both  A^^l  and  The  even  A  estimates  of 

and  odd  A  estimates  of  converge  quite  independ- 

ently  of  the  odd  A  estimates  of  and  th?  even  A  estimates 

of  3  It  ,  t  .  Boch  sets  converge  to  the  same  values,  so  it  is 

-  7  >  p  (O) 

necessary  to  compute  only  one  set.  The  set  which  begins  with 

is  chosen  and  so  as  to  reduce  the  number  of  coding  instructions,  per¬ 
sistence  is  chosen,  I  “  ^ 

The  convergence  test  may  be  satiefied  by  either  succeeding  values 
of  Bn  4  I  or  of  intervening  value  of  the  other  part 

is  then  also  accepted.  For  example,  if  13^4  |  within  tolerance 

of  ,  then  =.  • 

It  may  be  possible  to  design  iterative  procedures  for  solving 
Eq.  (74),  page  41,  which  converge  more  rapidly  than  does  Eq.  (76). 

In  this  respect  the  techniques  discussed  in  Section  10,  in  connection 
with  the  iterative  method  of  solving  the  finite- difference  analogue  of 

Poisson's  equation,  might  prove  successful.  However,  the  analysis  by 

10  r 

Young  does  not  apply  because  our  matrix  L  does  not  have  the 

necessary  property  (property  A). 
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7.  THE  COMPUTING  PROCEDURE  AND  OUTPUT 


According  to  Eq.  (42),  page  23,  and  the  marching  procedure,  the 
equations  which  must  be  programmed  are: 


where 


The  's,  'X  »  K  and  the  constants  required 

by  the  cvar  functions  are  precalculated  and  are  stored  in  the  computer 
A  consideration  of  the  types  of  profiles  to  be  investigated,  of 
computer  storage,  and  of  the  desirability  to  retain  any  symmetries, 
led  to  the  approximation  of  the  B  continuum  by  33  points,  32 
intervals. 

Regarding  the  units  of  dimension,  both  the  length  unit  and  time 
unit  shall  be  selected  from  the  straight  parallel  flow,  thus  keeping 


• 

r.*,  <  ■.'1  V  -/  C 


>  *-•  V  *.• 
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K/  variable.  Units  which  make 


U 

U' 


im/iL 


) 

5 


are  selected  for  each  profile.  In  atmospheric  considerations,  a  lower 
bound  for  the  time  unit  is  obtained  by  selecting  extreme  values  for  the 

A.  1 

shear.  In  the  horizontal  a  representative  high  value  is  10”  seconds”  , 
which  yields  a  time  unit  of  about  2  1/2  hours.  This  makes  the  length 
unit  corresponding  to  about  a  U  of  50  meters  seconds  450 
kilometers.  In  the  vertical  a  representative  high  value  of  the  shear 
is  10  ^  seconds  which  yields  a  time  unit  of  only  about  11/2  minutes. 
This  makes  the  length  unit  corresponding  to  a  U  of  50  meters 
seconds"  ,  4  1/2  kilometers. 

^  To  arrive  at  a  graphical  representation  of  the  evolution  we  use 

A(2,t)  +  i-BCz.t)  and  =  (X  Cz  ,i)  +  1/3 

where  A  ,  B  ,  OC  and  ^3  are  real  functions  of  Z  and  t.  . 

This  leads  to 


where 


■-  A(E,t) 

coti,  kx  - 

B(z„t)  kx 

=  Dte.i) 

COIL  [kx 

Ci.z.t) 

=  OC  (H,t) 

cotL.  kx- 

sx/ru  kx 

X(2,t) 

CertL  [kx 

-0u.t)]  , 

II 

o 

aA.ci<xAu  (7  B/A)  ^ 

X  = 

0  = 

QA.cjt'ayny  (•"^^(3^)  . 

*ln  a  frame  of  reference  moving  with  the  fluid  at  some  level. 


The  computer*  is  instructed  to  output  the  set  of  33  A's  and  the 
33  B's  for  each  time  unit,  that  is,  at  intervals  of  1/V  time  steps. 

A  second  code  is  later  used  to  compute  the  33  OC  's  and  the  33  /3  *8 
from  each  set.  According  to  Eq.  (33),  page  ZZ*  they  are  given  by 

K3(AJ,] 

m  J 

t  §  ^  *,(B)  . 

In  addition,  this  second  code  computes  a  figure  representative  of  the 
perturbation  kinetic  energy,  which  will  be  used  as  an  indication  of  the 
grov,'th  or  decay  of  this  quantity  with  time.  This  energy-  computation 
ic  based  on  a  linear  interpolation  of  between  points.  A  third 

code**  computes  D  and  (f)  from  each  A  ,  B  pair  and  X  and 
Q  from  each  OC  ,  /3  pair  according  to  Eq,  (9Z).  The  D  and  X 
retain  all  the  significant  digits  of  the  A,  B  and  OC  ,  /3  pairs, 
respectively.  The  ^  and  0  are  obtained  from  a  table-read 

_3 

routine  and  are  accurate  to  10  degrees. 


*  The  Standards  Western  Automatic  Computer  (SWAC),  Numerical 
Analysis  Research  (N.A.R. ),  University  of  California  at  Los 
Angeles. 

**  These  codes  are  on  file  at  the  N.A.R,  library. 


8.  APPLICATION  TO  AN  UNBOUNDED  HYPEKBOLIC- TANGENT 
PROFILE 

The  straight  parallel  flow  profile  is  given  by 

U  U)  =  -  tavnJv  Z* 

where  z  extends  to  ±  OO  . 

This  profile  is  one  which  is  acceptable  for  the  "contained"  treat  •’ 
ment  discussed  on  page  19.  By  choosing  a  sufficiently  deep  layer 
about  2=0,  U  (2)  can  be  made  as  small  as  desired  in  the 
neglected  exterior. 

This  profile  is  skew -symmetric,  that  is, 

U(H)  =  -U("2)  . 

I  I  j  1^  ® 

Hence  also  U  (z)  =  —  U  ("2)  ^  To  preserve  this  property  in 
our  space-discretized  system,  we  center  our  grid  at  the  level  Z  =  0 
with  16  intervals  above  and  below. 

It  is  usual  to  order  the  components  of  our  vector,  ^  ,  and  our 
set  of  equations  in  ascending  order  of  the  grid  point${that  is,  ascending 
M),  However  it  v/ill  serve  us  better  to  adopt  the  order 

M  =  -16,  -15,...  -1;  0;  16,  15,..,  1. 

Because  we  have  no  boundaries,  cL  of  Eq.  (45),  page  25,  is 

zero,  and  the  elements  of  C  >  which  we  shall  label  C..  ,  are 

~  M.j 

given  by 


C„,„=(tU„.u:K,) 

S  =(u:K3/a)K'”-^  jVM 


where  the  columns  are  labeled  as  the  rows,  in  the  order 
j  ::  -16,  -15,  ...  -1;  0;  16,  15,  ...  1. 

*  In  arbitrary  units,  U  (z)  =  U  tank  2  ,  and  we  are  free  to  add  an 
arbitrary  translation  to  our  frame  of  reference. 


By  virtue  of  Eq.  (95), 


c  .  =  0,  c  .  =  -c  .  . 

0,J  M.J 

To  exhaust  the  relationships  expressed  by  Eq.  (97),  we  write  the 
relationship  which  the  eigensolutions  must  satisfy  as  follows 


-lb - 1  0 


-16  / 


\-£ 


\  \  A 


which  is  the  same  as 


=  X 


= 


-p  -  %£  -  “V  " 

We  may  note  that  if  ^  ^  0  then  =  0  . 

It  is  apparent  that  one  of  the  latent  roots  is  =  0 

corresponding  vector  satisfies 

1^  p'  l%i\  K' 

\  ~  ~  7  \  / 

which  shows  this  vector  to  be  symmetric,  that  is. 
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Tht  *^ent  roots  in  the  case  of  skew -symmetric  profiles  occur  in 
±  pairs.  If  ^  ~  )  satisfies  £q.  (99),  then  also 

does  ~  ^  )  1  ^  5  ^  follows  from  Eq.  (99)  that  the 

satisfy 


sixteen  admissible  values  uf 


=  0 


(101) 


This  is  a  polynomial  in  with  real  coefficients.  Hence  an  ad¬ 
missible  is  either  real  or  is  one  of  a  complex  conjugate  pair. 

A  complex  conjugate  pair  of  X  's  give  rise  to  two  pairs  of  complex 
conjugate  roots,  X  =  yu  ±  iV  and  X=  “(yU.  ±  lv}  .  To 
each  such  pair  there  corresponds  a  complex  conjugate  pair  of  latent 
vectors. 

It  follows  from  the  preceding  paragraph  that  if  Eq.  (101)  admits 
a  real  X^  which  is  negative  (that  is,  X  =  ±  lV  where  V  is 
real),  then  the  corresponding  pair  of  solutions  may  be  expressed  by 


bV  ^  { P+  iQ,  0,  P-  iQ| 
-lv  ^  (P*  *-Q>  0)  P*^  iQj 


where 


and 


__  Q  are  real. 

The  profile  Eq.  (94)  has  already  been  studied  by  Garcia.  ^  He 
chose  as  the  principal  dependent  variable  and  investigated  the 

steady- state  differential  equation 


-Uci) 


^"(2) -k.'-4^(z)  +  U"(z)^(z)*0 


(102) 


(103) 


This  steady-state  condition  is  readily  obtained  from  Eq.  (8),  page  5, 
with  substitution  for  according  to  Eq.  (9). 

With  U  =  “  tank  z,  Eq.  (103)  has  the  solutions 

'1^,  -  Uu)  +  - 
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To  satisfy  the  boundary  conditions  (which  neither  one  of  the  above  can 
do)  Garcia  fitted  these  together  (about  the  level  0  )  arriving  at  the 
solution 

%  =  [t  »lU(H)l]  e (104) 

Since  ^  as  well  as  ^  must  be  continuous,  Garcia  found  ft  =  1 
to  yield  the  only  possible  stationary  wave,  no  matter  at  which  level 
the  two  solutions  are  fitted. 

For  ^  ^  Eq.  (104)  has  a  discontinuity  in  (that  is,  in 

■U.)  at  Z  -  0.  The  harmonic  field  which  is  compatible  with  this 
slipping  (sometimes  called  sliding  vorticity)  is 

=  Ck.- i/t)  (105) 


Removing  (jL  from  Eq.  (104)  Garcia  arrived  at  tlte  continuous  field 


l/li  +  |U(?) 


(106) 


which  is  no  longer  stationary.  From  the  initial  phase  speed  at  all 
levels,  he  found  that  for  k.  ^  i  an  amplifying  tilt  develops  and  for 
>  I  a  damping  ^ilt  develops,  This  indicates  the  imminent 
tendency.  Garcia  concluded  that  the  short  waves  are  stable  while  the 
long  waves  are  unstable  and  =  1. 

For  extremely  short  waves  (  h.  — ♦  oo  )»  the  determinant,  Eq.  (101), 
reduces  to  a  diagonal  with  elements  •  This  means 

that  the  components  of  ^  ,  that  is,  the  's,  are  each  advected 

by  the  straight  parallel  flow  at  their  level,  their  influence  on  other 
levels  being  negligible. 


♦According  to  the  Reynolds-Fjortoft  criterion.  See,  e.  g.  , 
Holmboe^  p.  11. 
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Thus  for  large  fl  ,  all  (k.)  are  real  and  positive.  In  the 
light  of  Garcia's  work  we  may  expect  one  of  the  16  values  of  (k.) 
to  become  sero  in  the  vicinity  of  k.  =  L.  This  value  of  X^  (K.)  may 
become  negative  as  k  decreases  further.  The  numerical  results  to 
follow  agree  with  this  conjecture. 

Garcia's  stream  function,  Eq.  (106),  is  taken  as  the  initial  field. 

It  is  compatible  with 

(^(2)=  -  E  (k  +  [tdLAliv,  *  j)  2  e  ,  (107) 

which  is  symmetric  and  real.  *  The  amplitude  factor  remains  free 
because  it  drops  out  of  the  homogeneous  linear  system. 

The  space  increment,  k  ,  was  chosen  to  be  0.2.  This  extends 
our  grid  from  z  *  -3.2  to  z  >  3.2.  Figure  1  shows  U(z)  and 
U  (z)  in  this  interval. 


Fig.  1.  The  hyperbolic  tangent  profile. 
*Such  an  L;iitial  value  may  prove  to  lack  generality. 
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Integrations  were  made  in  steps  X  =1/20  or  less,  from  t  =  0  ' 
to  t  *  15  for  the  values  k.  =  0.3,  0.5,  0.7,  and  2. 0.  The  last 
of  these,  k.  =  2.0,  was  apparently  stable.  The  others  lying  within 
0<  k-  <  1  were  unstable. 

In  each  of  the  three  unstable  cases,  a  single  X  eigensolution 
corresponding  to  X  =  IV  (  V  positive)  rapidly  dominated  the 
evolution.  The  growth  rates  were  determined  to  be 


k.  =  0.3  0.5  0.7 

V  =  0.  174  0.  184  0.  129 


(108) 


These  values  of  V  are  considered  to  be  accurate  to  within  2  percent.  * 
In  view  of  thia  wide  margin  the  correction  afforded  us  by  Eq.  (68), 
page  37,  is  trivial.  On  the  basis  of  these  values  the  solid  curve  in 
Fig.  2  was  drawn. 

The  exact  wave-number  at  which  the  negative  ^  becomes  zero 
is  rather  difficult  to  locate  by  integrations.  These  integrations  would 
have  to  be  extended  for  long  periods  to  detect  the  presence  or  absence 
of  slow  growth  rates. 


Fig.  2.  The  unstable  band. 


*More  accuracy  could  have  been  attained  by  carrying  each  integration 
further,  thereby  purifying  the  unstable  eigensolution. 
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The  broken  line  extension  of  our  curve  has  been  drawn  to  the  point 
k.  =  1,  V  =  0.  It  is  not  expected  to  be  much  in  error. 

It  seems  reasonable  to  assume  that  the  other  end  of  our  curve 
(approaching  infinite  wave  length)  will  asymptotically  approach  the 
dotted  curve.  This  dotted  curve  represents  the  instability  of  the 
kinked  profile  (the  dotted  curve  in  Fig.  1),  of  the  same  dimensions. 

3 

which  has  been  investigated  by  Holmboe. 

We  will  assume  that  except  for  the  one  value  of  which 

passes  through  zero  near  k.  =  1  and  which  is  negative  for 
0  <  k.  4  1  .  all  other  15  values  of  ^  Ck)  are  real  and  positive 
for  0  4  k.  4  oO  .  This  assumption  may  be  verified  by  numeri¬ 
cally  computing  all  admitted  by  Eq.  (101)  for  discrete  values 

of  k.  .  However  there  is  no  indication  why  the  assumption  should 
not  verify. 

Because  k.  =  0.5  is  about  the  most  unstable,  only  its  evolve- 
ment  from  the  initial  state,  Eq.  (107),  will  be  given  in  detail.  The 
evolutions  for  k.  =  0. 3  and  0. 7  are  basically  similar.  The 
evolution  for  k.  =2.0  showed  a  rapid  fall-off  in  energy  to  lees 
than  1  percent  of  the  initial  by  t  =  6,  at  which  time  gross-truncation 
error  became  evident. 

The  evolution  for  K.  =  0.  5,  wave  length  4“^^  ,  is  shown  in  the  fol¬ 
lowing  diagrams.  The  representation,  Eq.  (91),  is  used.  Figure  3 
shows  the  amplitude  and  phase  of  the  perturbation  vorticity  at  times 
t  =  0,  5,  10,  and  15.  Figure  4  shows  the  same  for  the  perturbation 
stream  function.  The  "energy"  growth  of  the  perturbation  is  shown  in 
Fig.  5. 

To  draw  the  total  streamlines  and  total  vorticity  isopleths  we 
must  prescribe  the  amplitude  of  the  perturbation.  We  may  prescribe 
it  arbitrarily  large  to  show  the  features  of  the  flow.  The  lateral  dis¬ 
placement  of  a  total  streamline  at  a  level  z  ,  consistent  with  the 
linearization,  is  given  by 

_  X  (z.t)  cenL  r  kx  -  0Cz,t)l 

U(i)  ^  J 


(109) 
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Fig.  3.  The  q -functions  D(z)  and  ^  (z).  long  wave  case 
(k  s  1/2)  at  times  t  «  0,  5,  10  and  15. 
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Fig.  4.  The  yr -functions  %(z)  A^nd  9  (s).  long  wave  case 

(k  s  1/2),  at  times  t  »  0,  5,  10  and  15. 
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Fig.  5.  Kinetic  energy  curve  for  K  =  1/2. 

♦ 

and  tha^  of  the  total  vorticity  isopleth  by 

_  .-Oig-ill  cc^  [lix-iicz.t)]  (110) 

These  exprsscious  apply  only  as  levels  where  the  denominator,  UCa) 

Uii 

(z)  ,  is  an  order  of  magnitiide  greater  than  the  numerator. 
Figure  6  shows  the  total  streamlines  (solid  curves)  and  total  vor¬ 
ticity  isopleths  (broken  curves),  at  selected  levels,  at  time  t  s  0 
above  and  t  =  15  below.  The  exaggerated  amplitude  of  the  per¬ 
turbation  has  been  reduced  by  a  factor  of  10  in  the  lower  diagram 
relative  to  that  of  the  upper  diagram. 

At  no  time  does  it  appear  as  if  the  vorticity  in  the  neglected 
exterior  has  become  significant. 

All  integrations  and  auxiliary  computations  carried  out  for  this 
skew -symmetric  profile,  with  symmetric  initial  values,  retained 
their  symmetries  to  the  last  digit.  Some  twenty  hours  of  high-speed 
computing  were  involved.  This  speaks  highly  of  SWAC.  * 

*See  footnote,  p.  51. 
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Fig.  6.  The  total  ctreamlinei  and  total  vorticlty  ieolines  of  selected 
selected  levels,  long  wave  case  at  t  3  0  above  and  t  ^  15  below. 
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9.  THE  SIGNIFICANCE  OF  THE  RESULTS 


We  have  crossed  three  bridges  to  get  to  our  numerical  results 
The  first  bridge  connects  the  governing  finite -amplitude  system  to 
the  linearized  system.  The  second  connects  the  linearized  system 
to  the  space- discretized  system.  And  the  third  connects  the  latter 
to  the  space-and-time  discretized  system.  If  we  are  to  consider 
that  our  numerical  evolutions  represent  first-order  {in  amplitude) 
approximations  to  the  evolutions  of  the  firdte-amplitude  systems 
then  we  must  reconcile  these  crossings. 

The  analysis  of  Section  5  has  shown  that  the  character  of  the 
system  is  unchanged  by  crossing  the  third  bridge  provided  we  use 
the  centered  fjrst-order  analogue.  The  eigensolution  ^  , 

where  A  =  yu  -r  I.  V  becomes  the  l  eigensolution 
The  vector  is  unchanged.  According  to  Eq.  (68),  X-f  l3  given  by 


+  tv- 


=.  —  OA-cto/P^ 

V 


l-R"- 


1 


1  +  +  VT 


1 f  -  'n 


where 

R'  =  Ct/if  (/  +  ■ 

It  should  be  noted  that  each  of  the  four  quadrants  of  the  complex  X 
is  altered  in  the  same  way.  That  is,  any  symmetries  in  the  occur¬ 
ence  of  eigenvalues  are  preserved. 

We  cannot  as  readily  or  as  completely  reconcile  the  second 
crossing. 


Let  us  consider  what  is  revealed  by  the  eigensolutions  of  the 
linearized  system.  The  eigensolution  £^(.Z,i)=  C^C*) 
corresponds  to  of  the  space-discretized  system.  Its 


introduction  into  Eq.  (8),  page  5 ,  leads  to 


■n 


for  the  unbounded  model.  It  has  been  imposed  that  be  real  and 

that  h.  be  real  and  positive. 

Equation  (112)  may  be  put  into  a  standard  form, 

(rA)qu)=  Ute)i«,f)-u%)  [  (jcr)d/- ,  <1 

-OO  ^ 


where 


I(2i)  ^  ^  Y? 

o,_»o+  o-  +U-2r 


The  proof  that  Eq.  (113)  is  a  valid  expression  of  £q.  (112)  re¬ 
quires  that 


oQ 

J  %{?)djr  = '^(2)  • 


For  sufficiently  small  positive  a,  I  (s,^)  «  0  except  for  a 
narrov/  region  of  half-width  2a  about  ^  Z  .  This  permits  us 
to  take  a  mean  value  of  *  for  this  region,  outside  of  the 

integral  in  Eq.  (115).  This  value  is  QCz)  •  It  follows  that 

r  T 

Icz,^)  =  c^(h)  Icz.Oci^  =  ^Cz)  • 


The  second  integral  may  be  found  in  tables.  -It  is  unity  for  all  Z 
and  positive  a.  It  should  be  noted  that  the  limit  is  not  taken  all 
the  way  to  a  =  0  but  rather  to  0-^  which  can  be  as  close  to 
zero  from  above  as  desired. 

The  form,  Eq.  (113),  classifies  our  integral  equation  as  a 
singular  homogeneous  Fredholm  equation  of  the  second  kind.  The 
determination  of  the  possessed  by  Eq.  (113)  for  the  profile. 


Eq.  (94),  for  the  real  positive  range  of  k.  is  beyond  the  scope  of  the 

present  work.  We  shall  deal  only  with  some  salient  points. 

We  may  begin  by  indicating  some  of  the  similarities  which  Eqs. 

(112)  or  (113)  has  with  the  space-discretized  system.  It  can  be 

demonstrated  directly  that  if  X  and  satisfy  Eq.  (112)  then 

also  do  and  «  where  the  subscript  o  denotes  the 

complex  conjugate.  Because  of  the  profile's  skew-symmetry  it  can 

also  be  demonstrated  that  if  X  =  X(z)  satisfy  Eq.  (112) 

then  also  do  —  ^  and  Q  =  XC-Z)  •  It  follows  that  the  two  eigen- 

^  T'Z. 

functions  corresponding  to  a  I  which  is  real  and  negative  (if 
such  is  admitted)  may  be  expressed  by 

r  PC2)+  ^  >  0 

LV,  (^(2)  =  \ 

^  P(-2)-LQt-s)  z  <  0 

f  P(2)  -  IQ(2)  W  E  >  0 

-  Iv  ,  q(2)  ==  i 

(  P(-2)+iQ(-E)  fir'u  z  c  0 

which  is  the  counterpart  of  Eq.  (102). 

The  coxinterpart  of  Eq.  (101)  is  also  readily  developed; 

OO  OO 

O  o 

where 


Yte)  =  ' 

=  ux- . 


5 


There  is  one  very  significant  difference  between  the  linearizec 
system  and  the  space-discretized  system.  It  has  been  shown*  that 
the  zeros  of  (T"  on  the  real  axis;  neutral  wave)  must 

coincide  with  the  zeros  of  U' (z)  .  For  the  profile,  Eq.  (94), 

there  is  only  one  level  where  U  (z)  is  zero,  and  that  is  at  the 
level  2  =  0  where  LJ  =  0.  Hence  the  only  value  of  T'  on  the 
real  axis  which  may  be  admissible  is  T  =  0.  Garcia  has  shown 
further  that  P  =  0  is  only  admissible  for  k.  =  1. 

With  this  knowledge  it  can  be  shown  that  the  eigensolutions  of 
Eq.  (112)  are  not  complete.  Equation  (112)  evaluated  at  Z  =  0 
gives;  “P  (^(o)  =  0  .  For  k.  1  ,  P  =  0  is  not  admissible. 

All  eigenfunctions  must  then  be  zero  at  Z  =  0  and  it  follows  that  w 
cannot  fit  a  value  of  (^(p)  ^  0  by  combining  eigensolutions  if 

k.  ^  1  . 

We  must  conclude  that  all  the  real  positive  X  ,  which  are  ad¬ 
mitted  by  Eq,  (101),  yield  space-discr etizsd  eigensolutions  which 
are  extraneous.  They  are  not  approximations  of  eigensolutions  of 
the  linearized  system  because  the  linearized  system  can  have  no 
such  counterpart. 

The  extraneous  eigensolutions  appear  because  the  space- 
discretized  system  must  yield  a  complete  set.  The  system  can 
evolve  all  arbitrary  initial  fields  with  the  aid  of  these  eigensolu¬ 
tions  but  the  approximations  necessarily  deteriorate  with  time. 

The  extraneous  eigensolutions  cannot  standalone;  either  they  rely 
on  gross-truncation-error  or  they  are  not  capable  of  satisfying  the 
boundary  conditions  in  the  limit  as  we  approach  the  continuous 
system. 

It  seems  reasonable  to  believe  that  the  real  negative  value  of 
Tv?”  (k.)  »  0*^  1  .  yields  eigensolutions  which  are  legitimate 

counterparts  of  corresponding  eigensolutions  of  the  linearized  sys¬ 
tem.  This  implies  that  Eq.  (117)  admits  only  one  value  of  P  OO 

♦See,  e.  g.  ,  Garcia^  p.  90,  or  Lin^  p.  118. 


which  is  real  and  negative  for  0  <  K  ^  ^  and  which  is  zero  and  termi¬ 
nates  at  K.  =  1  (Garcia's  solution);  and  that  there  are  no  F'  admitted 
for  k  >  1 

Havi  ig  accepted  the  two  eigensolutions  with  eigenvalues  ±  t.  V 
for  0  <  k  ^  1  as  true  eigensolutions  of  the  linearized  system,  it 
remains  to  determine  their  significance  in  the  finite -amplitude  system. 

It  is  well  known  that  the  total  ene^y  (kinetic)  is  separable  into  the 
energy  of  the  mean  flow  defined  by  Li  —  "T"  J'  IL  where  L_ 

is  the  wave  length  of  the  perturbation,  and  the  energy  of  the  remainder 
of  the  flow  (the  eddy).  This  separation  is  possible  because  the  space- 
integral  of  the  product  term  is  zero  by  definition  of  the  mean  flow. 


During  evolution  the  finite-amplitude  barotropic  mechanism  may  trans¬ 
fer  energy  between  the  mean  flow  and  the  eddy  but  the  sum  of  the  two 
partial  energies  remains  constant. 


The  linearized  system  violates  this  energy  conservation.  The 
perturbation  remains  sinusoidal  in  X  and  hence  remains  the  eddy; 
the  straight  parallel  flow  which  is  held  constant  by  choice,  remains 


the  mean  flow,  'I'hus  the  energy  of  the  mean  flow  is  constant  even  as 
the  eddy  grows. 


The  inclusion  of  the  nonlinear  part  of  the  mechanism  necessarily 
rights  this  situation.  The  straight  parallel  flow  may  still  be  maintained 
constant  by  choice  and  may  make  up  all  or  part  of  the  mean  flow.  The 
perturbation  does  not  then  retain  its  sinusoidal  dependency  although  it 
does  maintain  its  periodicity.  Consequently  the  evolving  perturbation 
is  made  up  of  the  eddy  and  an  addition  to  the  mean  flow. 

It  is  of  particular  interest  to  determine  how  the  mean  flow  evolves 
if  the  initial  profile  is  dynamically  unstable.  We  can  determine  this 


and  the  significance  of  the  unstable  eigensolution  by  making  a  finite- 
amplitude  integration.  As  initial  state,  we  take  the  straight  parallel 
flow  plus  a  self- excited  unstable  eigensolution  (logically,  the  most 
unstable)  with  a  small  but  finite  amplitude. 

We  shall  not  further  discuss  the  damped  eigensolutions. 


10,  THE  FINITE-AMPLITUDE  PROBLEM 
The  finite-amplitude  evolution  obeys 

where  =  V  ^  .  We  shall  continue  to  keep  the  straight  parallel 

flow,  U  =  -tank  z,  constant  by  choice. 

Only  one  integration  is  carried  out.  The  perturbation  is  periodic 
in  X  with  period  4Tl  .  The  initial  perturbation  is  based  on 

£^(X,Z)  =  D(2)  CfftL  ,  (120) 


where  k  |:  1/  2  and  D  (2)  and  (2)  are  given  in  Fig.  3  at 
t  =15  except  that  0  (H.)  is  smoothed  through  2=0. 

A  grid  is  introduced  having  272  independent  points  (16  columns  by 
17  rows).  The  17  rows,  at  intervals  of  TT/S,  extend  the  grid  from 
2  =  -Tf  to  'It  ,  inclusively.  The  columns  divide  the  period  4lT  into 

16  intervals  of  Tr/4  each. 

The  grid,  plus  an  ordering,  transforms  our  fields  into  vectors. 

The  actual  initial  is  given  in  Fig.  7  for  one  complete  wave  length. 

The  extremes  of  the  perturbation  vorticity  are  approximately  one-fourth 
as  big  as  those  of  the  profile. 

We  may  refer  to  a  grid  point  by  its  coordinates  relative  to  the 
central  leftmost  point  (0,  0): 


X.  =  ilT/4  , 

X  =  TaTT/8  , 


11  =  0,  1,  2  ...  15 

-m,  =  0,  ±1,  ...  ±8 


(121) 


Or  we  may  refer  to  a  point  simply  by  giving  its  address  (1,  m). 
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initial  q-iield  (amplified 


The  boundary  conditions  are  treated  essentially  in  the  same  way 
as  they  are  treated  in  the  unbounded  linear  integrations.  However, 
instead  of  using  a  Green's  function'integral  type  solution  over  the 
whole  grid,  we  use  it  only  to  determine  the  values  along  the  two  rows 
Tr\,  =  ±  8.  In  the  interior  we  may  then  use  the  more  rapid  iterative 
inversion  of  the  differential  operator  analogue. 

Any  particular  grid  point  (  )  together  with  the  corres¬ 

ponding  point  in  each  of  the  other  periods  constitutes  a  vortex  row 
with  uniform  strength,  5  ,  and  uniform  interval,  L  =  4Tf  .  The 
strength  of  each  vortex  is 


SA 


where  8P\  =  (  'n’/4)(  1^/8).  The  stream  function  which  is  compatible 
with  such  an  unbounded  and  infinitely  extended  single  row  is  given  at 
the  arbitrary  grid  point  JL,  yn,  by* 


A  in. 

4Tf 


JL 

\ze> 


(122) 


The  stream  function  at  the  point 
distribution,  is  then  given  by 


(  L  "TV  ),  due  to  the  entire 


TT(m-mo 

16  6 


(123) 


in  which  the  sum  is  taken  over  all  the  grid  points  as  given  by  Eq.  (121). 
It  may  be  noted  that  in  determining  the  stream  function  we  are  neglect¬ 
ing  the  vorticity  outside  our  region  of  attention,  as  we  did  in  the  linear 
integrations.  It  is  assumed  that  the  vorticity  on  the  outside  remains 
negligible.  This  is  something  which  will  have  to  be  checked  by  exam¬ 
ining  the  growth  at  the  extreme  rows. 


5 

♦See  Lamb  p.  224. 
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For  economy,  Eq.  (123)  shall  be  used  only  for  the  extreme  rows, 
yr<^  =  ±  8.  The  values  of  'P  ,  bo  determined  on  the  extreme  rows, 
then  play  the  role  of  boundary  values  for  the  interior  for  which  we 
shall  invert  a  finite -difference  analogue  of  .  At  an 

arbitrary  point  this  analogue  has  the  form 


(124) 


The  coefficient  of  has  intentionally  been  made  -I,  If 

1  =  16  then  i  +  1  is  taken  as  0,  and  if  j,  =  0  then  i  -  1 

is  taken  as  15  by  virtue  of  the  periodicity. 

The  interior  grid  is  given  some  order  A  =  A  [I,  r/v)  ■ 

Equation  (124),  applied  at  each  point  in  this  order,  then  yields  the 
system  of  simultaneous  equations 

C  =  d  .  (125) 

If  ’'i  ^  ±  7  then  d,>  =  0  *  ,  .If  m  =  ±  7  then  the 

known  boundary  value  which  appears  on  the  left-hiand  side  of  Eq.  (124) 
is  moved  to  the  right-hand  side,  and  ^ ^  ^  g 

Equation  (124)  is  a  partial  difference  equation  of  the  elliptic  type. 

It  is  convenient  to  invert  the  system,  Eq,  (125),  by  an  iterative  pro¬ 
cedure  known  as  the  method  of  successive  overrelaxation.  ♦  To  for¬ 
mally  exhibit  the  method,  we  first  express  C  as  the  sum  of  three 
matrices  ,x“d+L-i  .Iny,  nonzero  elements  appear 
only  above  the  main  diagonal.  In  l_  ,  nonzero  elements  appear  only 
below  the  main  diagonal.  The  matrix  I  is  the  identity  matrix. 

The  method  is  given  by 

^  '  +•  OJ  K  (126) 

*The  extrapolated  Liebmann  method  is  of  this  type. 
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where 


(12V, 

The  superscripts  in  parentheses  number  the  successive  approximations 
of  Vp  .  The  operation  is  carried  out  in  the  order  if  .  Hence  at  a 
point  ,  the  points  following,  r  >  ,  are  operated  on  by  U  i 

^  I 

and  the  points  preceding  which  are  once  later  estimates,  \  .  by  L  . 

If  we  substitute  ip  +  ,  where  vjy  satisfies 

Eq.  (125),  in  Eq.  (126)  we  find 

Tel‘>  =  T‘’'  ,  (128) 


where 


ax-  I  /  ' 

-  6oL  (I -to)  I  +  to  U 


(129) 


This  reveals  that  the  convergence  depends  on  the  eigenvalues  of 
which  are  given  by 


T 


T-xI 


=  0  . 


(130) 


The  convergence  rate  is  generally  dominated  by  the  eigenvalue  having 
the  largest  magnitude.  *  Let  this  magnitude  be  X  .  That  ,X  ^  1  is 
a  necessary  and  sufficient  condition  for  convergence.  The  rate  of 
convergence  may  be  expressed  by  —  X*\,  X  . 

We  are  still  free  to  select  the  ordering,  r  (  2.  )  ,  and  the 

value  of  CO  .  We  would  like  to  select  these  so  as  to  make  X  a 
minimum.  The  following  analysis  of  this  problem  is  based  on  a  de¬ 
tailed  paper  by  Young. 


♦This  has  been  discussed  in  Section  6  in  connection  with  Eq.  (79). 
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This  is  so  if  each  term  of  the  formal  expansion  of  the  determinant 
either  vanishes  because  one  or  more  of  its  elements  is  zero  or,  if  all 
its  elements  are  nonzero,  then  as  many  come  from  above  the  diagonal 
as  from  below.  The  above  equivalence  is  generally  latent  in  elliptic 
difference  equations  for  most,  but  not  all,  grid  selections.  The  equiva¬ 
lence  may  then  be  realized  by  choosing  one  of  a  more  or  less  large 
family  of  orderings  (called  "ronsistent"  orderings  by  Young). 

It  will  be  seen  that  all  consistent  orderings  lead  to  the  same  eigen¬ 
values  for  T  .  Young  then  proceeds  to  determine  the  best  co  for  this 
family  of  orderings.  The  possibility  remains,  however,  that  there  may 
exist  an  ordering  which  is  not  consistent  and  which  leads  to  a  smaller  A 
than  can  be  realized  with  a  consistent  ordering.  Young's  analysis  does 
not  rule  out  this  possibility. 

The  equivalence  is  latent  in  our  problem  if,  and  only  if,  we  divide 
the  X  cycle  into  an  even  number  of  intervals  greater  than  two — as 
we  have  done.  The  simplest  consistent  ordering  and  the  one  w'e  adopt 
is  then: 


r  =  8yn.  +  i/2  + 


57  for  i 

176  1/2  "  " 

56  1/2  "  " 

177  "  " 


even 

TVL 

odd 

odd 

If 

odd 

odd 

It 

even 

even 

n 

even 

(133) 


In  effect,  this  numbers  the  240  interior  points  as  one  would  count 
checkerboard  squares,  counting  the  light  squares  first  and  continuing 
with  the  dark  squares, beginning  each  traverse  at  the  lower  left. 
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We  may  now  proceed  from  Eq.  (132).  Let  us  denote  the  eigenvalues 

/^i  ’  *  *  *  ^L\C>  *  These 

depend  on  Eq.  (124)  and  the  grid  dimensions,  but  not  on  the  ordering. 

It  follows  from  Eq.  (132)  that  the  A  are  gi’/en  by 


An.+ 


(134) 


That  A  be  less  than  1  can  be  accomplished  for  some  t*)  only  if 
all  yW  satisfy  -1  <  Re:  yv  <  1.  If  in  addition  all  are  real, 

then  Young  has  shown  that  the  beet  ui  is  given  by 


'-^b  =  >  + 


I  4  (1-/)'/^ 


This  yields 


X  =  O).  -  I 


(135) 


(136) 


The  eigensolutions  cf  C  +  I  are  given  by 

vl;  (L,M)  =  aU.  ^  I  (m.48)  ^ 

where  L  and  M  take  on  the  values 
L  =  0,  1,  2,  ...  15 

M  =  -7,  -6,  ...  0,  ...  6,7. 


(137) 


The  corresponding  eigenvalues  are  given  by 


(L,M)  =  COV  ^  +  4  CCRL  (138) 
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Hence 


n  -  -I  i  ^  4  CcnL.-;^  (139) 

/  5  ^  to  _ 

-  0. 3846J 

For  this  value  of  yCl  ,  Eqs.  (135)  and  (136)  yield 

=■-  \.^S\ 2  A  -  0.(oJl2  (140) 

Determining  the  rate  of  convergence  by  —  l/i\  X  shows,  in  this 
particular  application,  that  using  the  best  Cj  accelerates  convergence 
by  a  factor  jf  12  over  and  above  the  factor  of  2  gained  by  relaxing  the 
grid  points  successively  with  Cu  =  i  in  a  consistent  order. 

The  initial  stream  function,  ^  ,  as  given  in  Fig.  8  was  deter¬ 
mined  in  the  manner  described,  from  the  initial  vorticity  given  in 
Fig.  7.  The  relaxation  w«  s  continued  until  Eq.  (124)  was  satisfied  at 
each  point  within  the  exacting  tolerance  of  2  2.  68  x  10  ,  in 

the  same  units. 

11.  THE  INTEGRATING  PROCEDURE  (NONLINEAR  PROBLEM) 

We  shall  take  up  the  problem  of  the  time-discretization  as  the 
first  stage  in  arriving  at  a  complete  numerical  analogue  for  our  non¬ 
linear  system,  Eq.  (119),  page  68.  This  will  be  follow  ed  by  the 
space  discretization  which  was  begun  in  Section  10  (for  (jj  ), 

The  system,  Eq.  (119),  is  of  the  first-order  in  time  and  may  be 
regarded  as  a  special  case  of  the  generalization 

^  X,  =  t  =  (Ml) 

where  is  one  of  a  finite  number  of  fields.  The  F"  's  at  any 

instant  are  space  dependent  only. 


We  shall  again  use  'C  to  denote  our  time  increment  and  mark  off 
the  instances  t  =  0,  TT ,  LX  ,  ...  (hJ  -  1) X  ,  MX  ,  (M  +  1)X  ... 

As  the  time-discretized  form  of  Eq.  (141)  we  offer 


where  the  bar  implies  the  mean  for  the  incremei.t,  that  is, 


(143 


Linearized  on  any  zero-order  stationary  state,  Eq.  (142)  reduces 
to  the  same  formula  used  in  our  linear  integrations.  For  its  solution 
we  shall  have  to  resort  to  a  method  of  successive  approximations 
(reiterations).  Because  of  the  nonlinearity  the  problem  of  convergence 
becomes  largely  a  matter  of  trial  and  error. 

For  propagated  effects  the  domain  of  dependence  (region  of  influ¬ 
ence)  using  Eq,  (142)  may  be  all  inclusive  for  any  X  even  after 
space -discretization. 

The  form,  Eq.  (119),  of  our  system  was  arrived  at  from  the 
governing  equation, 


^^Q=-VvQ, 

where  "V  =  J  >«■  V  $  and  Q  =  ,  by  explicit  refer¬ 

ence  to  the  zero-order  straight  parallel  flow  which  is  kept  constant 
by  choice; 


(14-j 


V=  Ui-tvy  Q  = 


(14f 


We  may  instead  refer  to  some  other  state  which  we  shall  call  the 
s  field  and  denote  by  subscript  s; 


V  =4  Q  =  Qs**  ^ 


(14 


which  satisfies 


X  -  ^  Qs=  0  • 


(] 


This  gives 

=  -  (\{  v)^-(yOs)-v-v  v<^  .  (1 

We  may  postulate  that  for  a  given  region  and  particular  time  interval 
there  exist  an  s  field  which  minimizes  in  some  prescribed  manner 
the  contribution  of  the  nonlinear  term  of  Eq.  (148).  The  amount  of 
linearity  exhibited  in  this  manner  may  be  called  the  inherent  linearity 
during  that  interval. 

The  evolution  of  the  large-scale  fields  of  motion  in  the  atmosphere 
probably  exhibit  a  high  degree  of  inlierent  linearity  even  up  to  periods 
of  a  day  or  so. 

The  degree  of  inherent  linearity  during  one  time-step  (interval) 
seems  pertinent  to  the  problem  of  convergence  of  the  reiterative 
process.,  It  may  be  tliat  if  90  percent  or  sc  of  the  change  comes  from 
the  linear  terms,  then  the  convergence  is  essentially  governed  by  these 
terms.  This  is  uncertain;  but  the  smallness  of  V  must  certainly  be 
significant.  It  may  be  noted  that  even  in  the  limit,  T  — ♦  0  ,  the 
mechanism  generally  remains  nonlinear. 

Another  significant  feature  of  the  implicit  formula,  Eq.  (142) 
follows:  Applied  to  the  barotropic  equation,  Eq.  (144),  this  formula¬ 
tion  preserves  the  conservation  of  total  kinetic  energy  which  is  inherent 
in  the  system  for  a  closed  region.  It  follows  that  the  time-truncation 
error  remains  bounded;  hence,  this  form  of  time -discretization  does 
not  introduce  any  computational  instability. 

The  conservation  of  total  kinetic  energy  is  derived  from  Eq.  (144), 
in  a  closed  region  ( vp  = the  boundary),  in  the  following  manner. 
Multiply  Eq.  (144)  by  and  integrate  over  the  region; 

V'  V  Q  dLo.  • 


To  show  that  the  right-hand  side  is  zero,  we  substitute 


i{)V-^0  =  V- (V0'}^)“O'P'7'V-QV-v^  , 

where  the  second  and  third  terms  are  zero  everywhere.  Then  by 
Gauss'  theorem 

A 

which  vanishes  because  there  is  no  outflow. 

This  leaves  us  with 

i  .  0 

A 

^  A 

A  .  < 

The  first  of  the  two  integrals  on  the  right-hand  side  is  zero  because 

A 

and  the  total  contained  vorticity  remains  constant.  Thus  we  have 
shown 

A 

The  proof  for  the  time-discretized  form 

S7^\b  -v^vl;  (149) 

+  »  rr9  ^ 
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is  similar.  We  multiply  by  Y  integrate  over  the  area.  As  before 

the  right-hand  side  drops  out  and  we  are  left  with 


A 

?  I K.;  +  S  ^  •  S'  r  ^^n)]  ^  =  0  • 

A  A 

The  second  integral  is  zero  because 

[f  (vt.r  f 

=  \  i ('"'  t.r  J-  =  tS'  ^ °  • 

A  A 

This  leaves  us  with 

\  da  =  \ 

A  A 

which  was  to  be  proven. 

The  preservation  of  this  inherent  property  after  space-discretiza¬ 
tion  leads  to  an  analogue  which  may  be  rather  difficult  to  use.  To 
achiev  .  it  we  multiply  Eq.  (149)  by  and  transform  the  equation 

into  the  form 

This  is  then  space-discretized  in  such  a  nnanner  that  the  summation 
of  the  term  on  the  right-hand  side  cancels  over  all  interior  grid  points. 
The  boundary  conditions  must  have  the  bracketed  term  zero  at  the 


boundary  grid  points. 

In  practice  this  analogue  must  then  be  solved  by  a  reiterative 
process.  This  may  prove  to  be  difficult. 

Because  of  the  nature  of  our  boundary  conditions,  we  have  not  fol¬ 
lowed  through  with  the  total  energy-conserving  space-discretization. 
Instead,  for  space  we  have  used  more  corventional  finite-differencing. 
We  have  already  shown  in  Section  10  how  and  the  boundary 

conditions  are  to  be  handled. 

Our  governing  equation,  Eq.  (119),  page  68,  time-discretized  has 
the  form 


(15i: 


This  is  .to  be  applied  to  our  grid  which  at  a  representative  interior 
point  has  the  spacing: 

H 

0 

K 

V/  o  C  ill  c  E 

K 

o 

5 


Compass  directions  are  used  to  label  neighbouring  points  relative  to 
the  arbitrary  point  at  which  Eq.  (151)  is  being  applied.  We  choose  to 
use  the  analogue 


i\  -  Vj  (/3  ^  j)  -  (.%-  J  -  %) 


where  OC  =  2  k  U|y^  .  /S  =  2  fv  and  \\  =  Tf/8. 

We  have  one  such  equation  for  each  interior  point.  The  vorticity 
changes  at  the  upper  and  lower  boundaries  may  be  obtained  by  a 
linear  extrapolation  of  the  changes  from  the  interior,  or  one  may 


(152 
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attempt  to  keep  the  values  of  constant  there  in  the  hope  that  both 

the  values  and  gradients  at  these  fringes  will  remain  small  anyway. 
We  will  try  the  latter  method  and  will  check  the  validity  of  this  pro- 
cedui-e  from  our  results. 

The  set  of  equations,  Eq.  (152),  for  all  interior  points  may  be 
packaged  into  the  representation 


This  analogue  must  also  be  solved  by  reiterating.  We  shall  attempt 
the  same  procedure  used  in  the  linear  integrations: 


where 


(») 


=  xr 


(M 

N+' 


,  (l+t 

U) 


is  determined  from 


by  the  method  of  Section  10. 


The  iterations  are  repeated  until  4.  =  S  determined  by  *3.^ 

t  ^  lv  'TKJ'i  \  T  N  ^ 


in  that  both  are  compatible  with  within  the  prescribed  toler¬ 

ance  at  all  interior  points. 

As  a  first  guess  for  ^  j  for  j  we  choose  linear 

extrapolation  in  time. 


+  N+1  Jn-1  +  i  Jn-I 


except  for  the  first  time-step  where  we  take 


We  choose  TT  =  1/  10.  It  is  hoped  that  the  inherent  linearity  will 
be  quite  large  for  so  small  time-steps  ^depending  as  it  does  on  the 
particular  evolution)  and  that  Eq.  (153)  will  converge  essentially  as  if 
it  were  linear. 

The  determination  of  Yk1-vI  from  +  >  requires  the  de¬ 
termination  of  on  the  upper  and  lower  boundaries  by  Green's 

method,  followed  by  the  successive  overrelaxation  of  the  interior 

C^) 

points  beginning  with  Tn  +  i  first  guess.  We  desire  compati¬ 
bility  of  +  I  fo  within  quite  small  tolerance,  but  the 

successive  estimates  of  these  vectors  need  not  be  so  highly  compati¬ 
ble.  It  would  seem  wasteful  to  make  them  so. 

In  an  attempt  to  accelerate  convergence  to  +  l  and  » 

the  following  scheme  is  incorporated  into  the  relaxation.  After  each 
compilete  grid  traverse  in  a  series  of  traverses,  the  tolerance  is 
raised  by  multiplication  with  a  prescribed  factor  greater  than  one. 
Hence,  the  more  traverses  it  takes  to  determine 

greater  its  incompatibility  with  tolerance  is  then  reset 

at  its  original  small  value  to  begin  the  next  series  of  traverses.  Thus, 
if  convergence  (all  within  tolerance)  is  found  on  one  traverse  (which 
signals  A~  =  b  )  then  I  and  +  \  have  been  found  to  the 

desired  prescribed  initial  tolerance.  This  scheme  considerably 
shortened  test  conrvputations  in  which  a  factor  of  1,25  was  used. 

12.  THE  RESULTS 

The  initial  perturbation  field  which  is  superimposed  on  the  un¬ 
bounded  hyperbolic -tangent  profile  is  given  in  Figs.  7  and  8,  pages 
69  and  75,  respectively.  Computations  have  been  carried  out  which 
take  the  evolution  to  t  «  5.  *  Outputs  were  obtained  every  five  time- 
steps,  that  is,  at  intervals  of  "t  =  1/2.  Figures  9  and  10  show  the 
fields  at  t  =  2  and  Figs.  11  and  12  at  t  »  5.  Six  more  decimal 
digits  were  carried  and  outputed  beyond  those  shown  in  the  figures. 

♦See  footnote  p.  51. 
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Fig.  10.  The  y- -field  at  t  «  2  (amplified  by  10  /I.  6). 
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Fig.  IZ.  The  \jr  -fielH  at  t  =  5  (amplified  by  10  /I.  6), 


It  is  notable  that  the  position  and  growth  develops  essentially  as 
predicted  by  the  linear  integrations.  This  shows  the  significance  of 
the  unstable  characteristic  yielded  by  the  linear  theory. 

What  we  call  the  perturbation  is  now  composed  of  an  eddy  and  a 
mean  (averaged  along  X  )  which  when  added  to  our  profile  gives  the 
mean  flow.  An  examination  of  Fig.  11  re^'eals  how  the  mean  flow  is 
being  modified.  The  centers  of  positive  perturbation  vorticity  have 
moved  inward,  toward  the  central  level,  and  the  centers  of  negative 
vorticity  have  moved  outward,  away  from  the  central  level. 

The  profile  itself  has  only  negative  vorticity  with  a  maximum  at 
the  central  level.  Thus,  the  eddy  is  altering  the  mean  flow  by  trans¬ 
porting  mean  negative  vorticity  away  from  the  central  level. 

The  perturbation  fields  at  t  =  5  are  averaged  along  x  to 
show  the  following  means. 


m 

8,  -8 

0 

4 

7,  -7 

-  22 

4 

6,  -6 

-  70 

-  1 

5,  -5 

-205 

-  15 

4,  -4 

-497 

-  62 

3,  -3 

-  593 

-186 

2,  -2 

-229 

-400 

1.  -1 

+860 

-650 

0 

+  1,  524 

-768 

These  values  of  and  have  been  amplified  by  10^/ 1.6. 

The  removal  of  these  means  from  the  perturbation  yields  the  eddy 
shown  in  Figs.  13  and  14.  Figure  15  shows  how  these  perturbation 
means  have  altered  the  mean  flow,  tending  apparently  to  stabilize  it. 
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and 


It  Is  quite  noticeable  that  the  error  due  to  .  ^ce-discretization  is 
becoming  more  serious  as  the  evolution  progresses.  This  is  one  reason 
why  the  computations  have  not  been  carried  further  with  our  present 
grid.  Another  reason  for  not  continuing  is  that  our  initial  value  lacks 
generality.  We  have  placed  the  constraint  of  periodicity  of  period  411 
on  the  evolution.  Asa  demonstration  of  what  may  develop  due  to 
dynamic  instability  and  the  presence  of  noise,  we  have  gone  far  enough. 

13.  CONCLUDING  REMARKS 

The  instability  mechanism  is  pertinent  to  a  number  of  meteorologi¬ 
cal  problems  on  all  scales.  However,  a  more  comprehensive  approach 
to  any  of  these  problems  generally  includes  the  effects  of  other  meteoro¬ 
logical  parameters. 

The  results  of  Section  12  demonstrate  the  modification  of  motion  on 
one  scale  by  the  formation  of  eddies  on  a  lower  scale,  the  scale  of  the 
eddies  being  dynamically  determined  by  the  causative  gradients.  It  has 
clearly  been  revealed  that  such  modification  cannot  be  handled  by  an 
empirical  eddy  viscosity.  The  profile  we  have  examined  contains  a 
vorticity  maximum  which  is  flattened  by  the  eddies.  If  the  profile  had 
a  vorticity  minimum  instead,  the  rate  of  smoothing  would  not  at  all  have 
progressed  so  rapidly  because  such  a  profile  is  dynamically  stable  and 
thus  does  not  favor  eddy  formation. 

Consider  the  large-scale  vertical  structure  of  the  atmosphere.  At 
a  single  locality  we  may  ignore  horizontal  variations  at  all  levels  if  it 
can  later  be  shown  that  these  variations  are  negligible  on  the  scale  of 
the  developing  eddies.  A  comprehensive  approach  to  the  stability 
problem  might  then  be  developed  which  takes  account  of  variations  with 
height  in  the  large-scale  density,  horizontal  velocity,  and  humidity 
distributions.  It  might  also  be  wise  to  include  the  modifying  effect  of 
the  local  variation  with  height  in  the  large-scale  fields  of  vertical 
motion  and  horizontal  divergence. 

The  eddies  which  form  due  to  the  vertical  structure  of  the  atmos¬ 
phere  are  not  only  significant  because  they  are  associated  with 
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convective  clouds  and  turbulence  but  also  because  they  modify  the  vertical 
structure. 

On  the  whole,  the  lapse  rate  must  be  the  most  significant  parameter 
in  determining  vertical  dynamic  stability  of  the  atmosphere.  Even  though 
it  can  be  shown  by  a  proper  selection  of  units  that  the  eddy  in  our  study 
grew  from  having  a  maximum  vertical  component  of  about  four  meters 
per  second  to  twelve  meters  per  second  in  less  than  eight  minutes,  this 
in  itself  is  not  significant.  Our  thinking  must  be  modified  by  other  con¬ 
siderations. 

We  must  also  consider  the  rate  at  which  the  dynamic  instability  of 

the  vertical  structure  is  being  built  up  by  large-scale  processes  in¬ 
cluding  surface  heating.  As  the  vertical  structure  crosses  thresholds 
of  instability  how  much  energy  becomes  available  to  the  eddies  ?  And 
are  the  eddies  capable  of  modifying  the  vertical  structure  rapidly  enough 
to  keep  the  instability  at  the  threshold?  What  then  is  the  strength  of  the 
eddies?  These  are  questions  which  may  be  answerable  by  numerical 
integration  studies. 

Knowledge  of  how  vertical  structures  are  modified  by  eddies  con¬ 
stitutes  a  prerequisite  for  making  dynamical  numerical  weather  predic¬ 
tions.  Consider  a  model  which  discretizes  the  vertical  structure  by  a 
number  of  levels,  and  which  has  grid  spacings  of  the  order  of  hundreds 
of  kilometers  in  the  horizontal.  If  no  allowance  is  made  for  the  vertical 
adjustments  by  eddies  of  a  smaller  scale  than  the  grid,  ridiculous  pro¬ 
files  and  lapse  rates  may  be  developed.  However,  if  the  mechanism  is 
understood  and  stratified  empirical  techniques  have  been  developed, 
then  not  only  could  the  vertical  structure  be  adjusted  to  keep  the  evolu¬ 
tion  on  the  right  track,  but  alio  useful  predictions  of  turbulence  and 
convective  cloud  formations  may  result. 

It  is  hoped  that  such  investigations  will  be  undertaken.  The  numer¬ 
ical  techniques  to  be  used  are  themselves  in  need  of  further  develop¬ 
ment. 
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If  the  behavioral  properties  of  a  differential  system  are  not  known 
and  are  being  numermally  investigated,  then  finite -difference  approxi¬ 
mations  and  the  order  of  truncation  error  are  of  primary  interest,  and 
the  concept  of  the  limit  is  essential.  However,  in  practice,  in  numeri¬ 
cally  integrating  partial  differential  systems  such  as  those  which  model 
atmospheric  circulations,  the  precept  that  we  are  bound  to  finite  dif¬ 
ferences  generally  exists.  The  differences  must  be  the  larger  the 
faster  we  wish  the  integration  to  proceed  relative  to  real  time,  and  the 
concept  of  the  limit  as  our  increments  tend  to  zero  does  not  enter. 

In  practice  we  are  concerned  with  the  complete  system  of  finite- 
difference  equations  as  a  numerical  analogue  of  the  complete  differen¬ 
tial  system. 

The  procedure  is  to  construct  and  modify  a  numerical  analogue, 
having  finite  increments,  until  its  behavioral  properties  resemble  as 
closely  as  possible  those  of  the  given  differential  system.  This  re¬ 
quires  a  pretty  good  understanding  of  the  differential  system. 

For  some  systems,  exact  numerical  analogues  can  be  derived. 

To  illustrate  this,  let  us  take  another  look  at  the  linear  system 

i  ^(t)  =  -  c  C  Q  Ci)  (156) 


We  arrived,  in  Section  5,  by  what  we  shall  call  the  finite-difference 
approach,  at  the  analogue 


(157) 


which  may  be  written  in  explicit  form. 


q  =  D  q  =  D  q  , 

TN-*-!  ^  TM  ~  TO 


(158) 
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where 


D  =  [lU  iTC)'(zl-itC  . 


(159) 


Our  analysis  showed  us  that  this  analogue  has  similar  properties  to 
the  differential  system  but  that  the  prediction  deteriorates  with  length 
of  forecast. 

If  we  abandon  for  the  moment  the  finite-difference  approach,  we 
can,  by  a  full  appreciation  of  the  properties  of  the  differential  system 
as  revealed  by  its  eigensolutions,  arrive  at  a  numerical  analogue 
which  exactly  expresses  the  evolution.  In  terms  of  the  eigenvectors 
and  eigenvalues  of  the  matrix  L/  • 


(160) 


gives  the  evolution  of  the  initial  vector 


intervals  “X  ,  If 


we  now 


construct  a  matrix  M  which  has  these 


same  eigenvectors  but  has  eigenvalues  ,  then 


q  »  M  q  =  M  q 

I’M  to  css  ts-l 


(161) 


which  expresses  £q.  (160)  exactly,  no  matter  the  magnitude  of  T 


*It  may  be  noticed  that  the  expression,  Eq.  (16 1),  resembles  a  linear 
regression  formula  in  which  the  same  set  of  parameters  enters  both 
as  predictor  and  predictand.  If  it  is  justifiable  to  assume  that  the 
set  of  parameters  is  governed  by  an  expression  such  as  Eq.  (156) 
then  from  an  analysis  of  the  matrix,  wUch  lias  been  determined 
statistically,  one  may  construct  both  the  tendency  matrix  and  a 
matrix  to  be  used  for  any  other  time  interval.  Some  difficulty,  how¬ 
ever,  will  be  encountered  in  interpreting  those  stable  eigensolutions 
which  have  undergone  more  than  one  oscillation  in  the  data  interval. 


We  have  demonstrated  that  for  certain  simple  differential  systems 
exact  (except  foi'  random  round-off)  numerical  analogues  are  possible^ 

We  have  also  indicated  that  we  should  not  be  constrained  by  the  finite- 
difference  approach  and  the  limit  concept. 

When  dealing  with  nonlinear  partial  differential  systems  the  problem 
becomes  rather  nebulous.  We  have  no  eigensolutiona  to  guide  us  and  we 
may  know  very  little  about  the  s>stem.  An  attempt  should  be  made  to 
learn  as  much  as  possible  about  the  behavioral  properties  of  the  systems. 
If  limit  cycles  are  understood  these  could  perhaps  be  modeled  by  the 
analogue. 

With  little  else  to  guide  us  in  complicated  systems,  then  perhaps 
the  best  we  can  do  is  to  attempt  to  preserve  physical  continuity  princi¬ 
pals  which  may  be  inherent  in  the  system.  An  example  of  this  is  £q. 
(150),  page  80,  which  attempts  to  preserve  the  conservation  of  total 
energy  in  a  closed  system.  It  is  not  yet  clear  how  many  such  princi¬ 
pals  one  can  preserve  simultaneously. 

This  last  approach  may  lead  to  conaidei'abls  success  for  General 
Circulation  models  (models  which  contain  atmospheric  forcing  and 
friction  terms  and  which  are  integrated  to  times  far  removed  from  the 
initial  values).  The  problems  are  many. 
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